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TERMINAL DESCENT SIMULATION STUDY 

By Richard F. Broderick, Michael K. Mann, 
Richard D. Moog, Billie W. Preston 
Parker S. Stafford, and Henry C. Von Struve 


SUMMARY 

This report documents the development of a Martian soft lander 
terminal descent simulation program. This program is designed 
for use relative to the parachute and terminal descent phases of 
the Mars Lander (soft landing) mission; however, it can also be 
used for other phases with slight modifications. The parachute- 
phase simulation was developed under this contract, but the 
terminal-descent-phase simulation was developed by Martin Marietta 
Corporation independently. Instructions for using the program are 
provided in a separate user's guide (ref. 1), which was also pre- 
pared under this contract. At the same time the simulation pro- 
gram was being developed, studies were conducted on several land- 
ing radar mechanizations to evaluate their performance character- 
istics in connection with the soft lander mission. These studies 
were performed using the terminal descent program and Monte Carlo 
simulations of two types, static and time-correlated. 


I . INTRODUCTION 

The delivery of a scientific payload to the surface of the 
planet Mars, with a soft-lander capsule, is a complex mission 
that includes an aerodynamic -entry phase, a parachute phase, and 
a terminal-descent, retro-engine phase. The terminal descent 
system must accomplish navigation, guidance, and control functions 
to position the capsule just above the surface with the prescribed 
velocity and vehicle attitudes. These functions must be performed 
regardless of the atmospheric winds, atmospheric densities, surface 
slopes, or other environmental variations encountered during the 
descent . 

Lander system studies have shown that multiple, differentially- 
throttled, monopropellant engines are the most efficient retro- 
engines and should be used for propulsion during the descent. 



For effective navigation, the relative range and velocity must 
be measured by Lander sensors during the terminal descent. These 
sensors are assumed to consist of a multiple-channel Doppler ve- 
locity radar and a radar altimeter that measures the distance 
to the surface. But the composite effect of the multiple-stage 
descent, the environmental variations, and the mechanization of 
the terminal descent system establishes a non-analytic performance 
problem. This creates the need for a six-degree-of -freedom (6D0F) 
simulation program to be used for parametric performance investi- 
gations. Since the initial terminal descent conditions are im- 
portant parameters, the 6D0F program should be capable of simu- 
lating the parachute phase as well as the terminal descent phase. 

Such a program has been developed by modifying a computer 
program that was originally prepared for the Air Force Systems 
Engineering Group by Litton Systems, Inc. (ref. 2). The resulting 
program is a powerful tool for soft-lander performance evaluations. 
Part of the modification of the old 6D0F program was accomplished 
under this contract. Since the performance of the doppler ve- 
locity radar is critical to the success of the mission, the 6D0F 
program and Monte Carlo simulations were used to investigate 
several candidate radar mechanizations, to demonstrate the analyti- 
cal techniques, and to provide background data on radar perform- 
ance. 

The report is arranged to present a description of the proposed 
mission, a description of the terminal descent guidance system, 
the results of the radar study, and a general description of the 
6D0F program. In this way, the reader gains sequential knowledge 
of the need, usage, and operation of the program. Detailed user's 
instructions may be found in reference 1. 


II. MISSION DESCRIPTION 


Continuous ly- active Lander Capsule propulsion, guidance, and 
control (PG&C) functions are required from Orbiter separation to 
touchdown on the Mars surface (figure 1). One navigation, two 
guidance, and several attitude-control modes of operation are re- 
quired to implement the sequential functions and performance re- 
quirements that are shown in the figure. The proper operation of 
the system is critical to the success of the soft landing mission* 
therefore, it must be reliable and flexible enough to compensate 
for variations in atmospheric density, horizontal winds, surface 
elevation, and surface slopes in the landing area. 

A block diagram of the PG&C system is shown in figure 2. All 
control computations and descent sequencing signals are provided 
by the general-purpose digital computer. 

Functional operation of the system begins 15 hr before the 
Lander separates from the Orbiter. Ground commands energize the 
system, initiate a checkout mode, and update nominal mission con- 
trol parameters. These parameters include the deorbit impulse and 
attitude, entry attitude, and sequence signal times. During the 
preseparation period, the GCC checks the system and relays system- 
status information to the Space Flight Operations Facility (SFOF) . 
An attitude reference is established in the GCC while the Orbiter” 
is locked on the Sun and Canopus. 

Separation is initiated by an Orbiter timer. An attitude-hold 
mode is used while the Lander coasts away from the Orbiter. Dur- 
ing this 30-minute period, the Lander is oriented to the deorbit 
attitude. Then the deorbit engine is ignited with a timed signal 
and is shut down when the VRU measures the required impulse. 

After the deorbit thrust phase, the Lander coasts in an atti- 
tude-hold mode until 30 minutes before entry. At that time the 
Lander is oriented to the programed entry attitude. The HARA is 
energized before entry with a timed signal. Entry is sensed by 
the VRU at 0. 1-g deceleration, and the control mode is changed to 
rate damping in all three axes. When maximum deceleration is 
sensed, local vertical inertial navigation equations are initial- 
ized in the GCC on the bases of a priori knowledge of velocity and 
attitude and the altitude determined by the HARA. At a preset 
altitude of 23 500 ft, as sensed by the HARA, the Ballute is de- 
ployed and control is switched to roll rate damping only. The 
parachuse is deployed at 14 500 ft. Six sec later the aeroshell 
is jettisoned and the TDLR and LARA are energized. 
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Attitude hold 
Orient for entry 
Jettison deorbit 
Module 


Control deorbit impulse 
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Parachute descent 
Roll rate damping 
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Figure 1.- Major Steps Between Orbiter Separation and Touchdown 



Figure 2.- Block Diagram of the PG&C System 
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During the parachute phase, the local vertical inertial navi- 
gation estimates of altitude and velocity are updated in the com- 
puter with radar data. At a preset altitude of approximately 
3900 ft, as sensed by the widebeam LARA, the terminal engines are 
ignited and the parachute is jettisoned. For the rest of the de- 
scent, pitch and yaw are controlled by differentially throttling 
the three terminal engines, and roll is controlled with six jets 
on the Lander. During a 3-sec tipup maneuver after parachute re- 
lease, the Lander is rotated to align the thrust vector with the 
velocity vector sensed by the TDLR. The Lander drops at throttle 
settings of 25% until the axial control curve conditions are met. 
At that time, the engines are throttled up to 90% and the Lander 
performs a gravity turn descent aimed at 10 fps at a 60-ft alti- 
tude. A constant- velocity control mode is then invoked to an al- 
titude of 10 ft, at which time the engines are shut down and the 
Lander drops to the surface. 



III. TERMINAL DESCENT GUIDANCE AND CONTROL 


A. Guidance and Control Laws 

The guidance and control laws considered in this study repre- 
sent the result of a series of studies documented in references 3 
thru 6. A summary of the descent guidance laws follows. 

To^ accomplish a soft landing, the three components of the 
Lander s velocity measured with respect to the surface must be 
controlled as a function of altitude measured with respect to the 
surface . Gravity-turn steering has been selected as the means of 
controlling the lateral velocity components. This law is mechan- 
lze y rotating the vehicle about its pitch and yaw axes until 
the lateral body-axis velocities are zero. This causes the thrust 
axis to point along the total velocity vector. This steering law 
is simple to mechanize because local vertical sensing is not re- 
quired. Gravity causes the thrust axis to rotate toward the ver- 
tical as the velocity is reduced. An arbitrary roll orientation 
is maintained by using an attitude-hold mode during the descent 
If the components of vehicle velocity are denoted as u, v, and 
w along the roll, pitch, and yaw axes, respectively, then the 
steering signals are 


Pitch: a = G w/u 

c a 

Yaw: 6 = G v/u 

c p 

Roll : cp = 0 

c 

where G^ and G^ are the pitch and yaw gains, respectively. 

, T J e axial exponent of velocity, u, is controlled by modulat- 
ing the vehicle thrust to follow the preprogramed desired velocity/ 
altitude contour shown in figure 3. After a 2-sec engine warmup, 
the parachute is jettisoned and a 3 -sec maneuver phase is alloted 
for aligning the thrust and velocity vectors. The high-thrust 
phase is designed to command up to 90% of the available thrust; 
this will conserve propellant by minimizing gravity losses. The 
remaining 10% of available thrust is used for steering. The con- 
stant-velocity phase allows the vehicle control system to reduce 
the effects of dynamic control lags and errors that result from 

following the high-thrust phase contour before shutting down the 
engines . 


Terminal descent contour Propellant used, lb 



Velocity (fps) 

Figure 3.- Terminal Descent Contour 

The contour is designed for the highest velocity condition, 
includes winds resulting from the parachute phase, and is shaped 
to account for propellant -tank pressure blowdown and specific- 
impulse changes. Simulations have shown that this approach allows 
margin for any 3a combination of atmospheres, winds, and surface 
slopes (ref. 4). The minimum velocity resulting from the para- 
chute phase occurs in the maximum-density atmosphere with zero 
wind. The resulting trajectory shown has a long drop phase before 
intersecting the design descent contour for the high-thrust phase. 
The propellant required for this condition is less than that re- 
quired for the maximum-wind case. 

* The propulsion-system design parameters of blowdown ratio, 
maximum thrust, throttle ratio, and propellant loading are defined 
through descent simulations using the guidance laws described 
above. The thrust and blowdown ratio can be chosen to minimize 
the weight of the system (ref. 6). Once the maximum thrust has 
been established, the required throttle ratio is defined by the 
requirement to reach a Mars thrust -to-weight ratio of less than 
one during the constant -velocity descent phase. Thus 

(1/R 4- 0.1) F/W < 1 

where R is the throttle ratio, F is the total maximum thrust, 
and W is the weight of the Lander on Mars. The 0.1 term in 
equation (4) is the margin reserved for attitude control. 



The use of gravity-turn steering with a high thrust -to-weight 
ratio can result in high vehicle turning rates near the surface. 
These potential effects are minimized by using a constant -velocity 
descent phase and a blowdown propulsion system. 

Figure 4 shows a block diagram of the Terminal Descent Guidance 
and Control System. Several choices exist for mechanizing the sys- 
tem; these will be discussed later in this report. 



Note: a - Velocity angle of attack - pitch axis 
0 - Velocity angle of attack - yaw axis 
V £ - Velocity error 
G - Gain 


Figure 4.- Block Diagram of the Terminal Descent Guidance 
and Control System 

As shown in the figure, the vehicle dynamic motion is sensed 
by the TDLR, the altimeter, and the three-axis strapdown gyro and 
accelerometer systems. The data from these sensors are used to 
generate outputs for u, v, w, and H from the radar-aided 
inertial navigator. The navigator equations are given below. 

Velocity equations: 

G = A x - qw + rv + gA 13 + K u (u R - u) ( 5 ) 

v = A y + pw - ru + gA 23 + K v (v R - v) ( 6 ) 

^ = A z ~ pv + qu + gA^ + K w (w R - w) ( 7 ) 
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Altitude equations: 

Z = A 13 u + A 23 v + A 33 w 




H = H - Z 
o a 


H - H, 


In these equations, u, v, w are the body-axis velocity compo- 
nents; u„, v__, w are the velocity components determined by 

the radar system; p, q, r are the attitude rates; A^, A , 

A are the acceleration components; A 13 , A^, A 3 3 are the 

z 

computed direction cosines; H is the computed altitude; H q is 
the reference altitude; H is the altitude determined by the 

3 

altimeter; g is the Mars gravitational constant; and K^, K^, 


K , K 
v w 


are adjustable gains. 


The attitude-rate data used in these equations are obtained by 
differencing the consecutive gyro-pulse counts in the digital com- 
puter. The acceleration data are obtained similarly from the ac- 
celerometer pulses. The gains K^, K^, K^, and are keyed 

to the radar-tracker lock status indicators. When a particular 
radar data input |e.g., U R j is invalid, the corresponding gain 

j is set to zero. This allows the computation to proceed using 

inertial data. When the gains are not zero, the inertial estimates 
will track the radar data with a time lag of 1/K sec /K = K^, 

V etc )* 

Three-axis accelerometer data are required to implement these 
equations. During the parachute phase, the inertial navigator 
updates the estimates of u, v, w, and H, and the lateral- 
axis accelerometers measure lateral accelerations caused by wind 
gusts . 

The descent contour block shown in figure 4 is a curve-fit 
function of the velocity /altitude contour shown in figure 3. Be- 
cause the form of this function for flight software has not yet 
been defined, the M0D6MV program uses a table look-up routine to 
generate the desired velocity from H. This desired velocity is 
then compared with u to generate the velocity error, V R . The 
command signal to the inner acceleration control loop is 


+ G VI ‘ S 
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where is the acceleration command proportional gain, is 

the acceleration command integral gain, and S is the Laplace 
transform variable. 


This signal is compared to the acceleration sensed by 

the roll-axis accelerometer (assuming G^ = lj, and is then fil- 
tered and limited. The filter removes accelerometer noise and 
structural feedback; the limiter ensures that the axial control 
channel does not saturate the engine valves and prevent differen- 
tial throttling for attitude control. 

In the steering channels, a rate command inner loop ensures 
that the initial pitch and yaw maneuver can be accomplished at a 
controlled rate after the parachute is released. The command 
rate signal is limited so neither the gyro-torquing capability 
nor the radar-tracking-rate capability will be exceeded. Pro- 
portional plus integral feedback is used to allow complete nulling 
of a c and in the steady-state condition, even in the pres- 

ence of misalignment torques. The M0D6MV simulations have shown 
that these integrators should not be active until after the vehicle 
has accomplished the initial maneuver. Unless this is done, the 
integrated signal tends to counteract the proportional signal, 
and the maneuver response is poor. After the maneuver is over, 
the integrators can be started to provide the desired steady-state 
nulling. 


Stability analyses (ref. 6) have shown that lead-lag compen- 
sation is desirable in each of the inner control loops (pitch rate, 
yaw rate, and acceleration) to improve their stability and response 
This lead-lag compensation is accomplished by using a filter trans- 
fer function. 


k l ( s + w „) 


hi/ 5 -* s + w 


( 12 ) 


on each of the input lines to the engine 
figure 4. In equation (12), is the 

filter lead break frequency, and is 

quency. In the present design (ref. 6), 
300 rad/sec, and K = 6. 


command mixing block of 
filter gain, is the 

the filter lag break fre- 

W = 50 rad/sec, W = 

° P 
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The engine command mixing block consists of the following 


equations for a three-engine 

vehicle 



B i = G P i • 

Pe + Te 


(13) 

B 2 = -G p2 

• Pe + G 

y 

• Ye + Te 

(14) 

B 3 = ” G p2 

• Pe - G 

y 

■ Ye + Te 

05) 


where B^, , B^ are the mixed command inputs to the engine 

valves, Pe is the pitch-channel error signal, Ye is the yaw- 

channel error signal, Te is the acceleration-loop error signal, 

and G l9 G 0 , G are the mixing gains in the three engines, 
pi p2 y 

These equations correspond to the engine arrangement shown 
in figure 5. For a symmetrical arrangment 


Pitch 


Yaw axis 



Figure 5.- Engine Arrangment 


The moment arms are related by 


Ls =V|"l 1 


(16) 

(17) 
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( 18 ) 

(19) 


The moment and thrust equations are then 

M p " 

= |f 2 • L 3 J - |f 3 • L 3 j 

F = -F x - F 2 - F 3 (20) 

where M is the pitch moment, M is the yaw moment, F is 
P y 

the total thrust of the engines and F x , F 2 , F 3 are the indi- 
vidual engine thrusts. 


G pl’ G p2 ’ and 
coupling between the 


G 0 should be chosen to eliminate cross- 
p3 

control channels. We can do this by letting 


F. = B. 

l i 


( 21 ) 


and using equations (13) thru (13) and (18 thru (21) to construct 
the matrix equation shown below. 



Since by equation (16) 



we see that choosing 


G 


p2 


(23) 


will yield a diagonal matrix which decouples the control channels. 
Furthermore, to make the pitch and yaw channels have equal gain 
through the mixing and engine elements, let 


G 

y 


2 G P 1 


(24) 
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The dynamics of all the components in the terminal descent con- 
trol system are well known, with the possible exception of the en- 
gine-throttling dynamics. Preliminary engine-throttling tests per- 
formed by Martin Marietta with nitrogen-saturated propellant show 
that the engine equivalent- throttling- time constant is on the or- 
der of 15 msec or less. The recommended control loop mechaniza- 
tion can be made stable and responsive with any postulated engine 
dynamics (up to a time constant of 100 msec) because all the other 
components in the system have fast responses, including the avail- 
able throttle valves (fig 6) . 



Figure 6.- Frequency Response to Input Commands 

The engines operate in a blowdown mode, which means their 
thrust (and therefore the gain in the inner loops) decreases as 
a function of the propellant used. The control loops are designed 
to have sufficient stability margin and to be stable and respon- 
sive at all operating gains, without requiring gain changing as 
a function of the propellant-tank pressure. The radar-aided in- 
ertial control loop mechanization enhances the stability and the 
speed of response of the control loops because the main control 
signals do not have to be filtered through a low-frequency noise 
filter (see ref. 5). Figure 7, which shows the root loci for the 
inner and outer control loops in the radar-aided inertial mode, 
illustrates that all dominant closed-loop poles are stable and 
well-damped. This was confirmed in the 6-D digital computer sim- 
ulations performed on the control system. As evidenced in the 
plots of pitch rate and position versus time, the response of the 
attitude channels is fast and stable. The axial channel also* 
has a fast response. In these analyses, the valve transfer func- 
tion was defined by a second-order linear model, and the engine, 
by a first-order model. 
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(a) Pitch-Yaw Control Inner Rate Loop Root Locus 
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(b) Pitch-Yaw Control Outer Radar-Aided Inertial Loop Root Locus, 
High-Gain Case 


Figure 7.- Root Locus Plots for the Inner and Outer Control Loops, 
System Operating in the Radar-Aided Inertial Mode 
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(c) Velocity Control Inner Loop Root Locus 



(d) Velocity Control Outer Radar-Aided Inertial Loop Root Locus, High 
Gain Case 


Figure 7.- Concluded 
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B. Alternative System Mechanization 


The block diagram in figure 4 showed a system that employed 
the outputs of a radar-aided inertial navigator for steering and 
descent-contour guidance. Alternatively, the outputs of the ra- 
dar |u R , v R , and w R j and the altimeter could be used 

directly for this purpose. This direct method is called opera- 
tion in the radar prime mode. Stability analyses (refs. 3 thru 6) 
have shown that for operation in the radar prime mode, special 
control-loop compensation is required to counteract the effects 
of having the radar noise-smoothing filters directly in the steer- 
ing channels. The compensation subsystem, which uses integrals 
of the gyro rate outputs q and r, is shown in figure 8. t 
is the time constant for the filter. This subsystem is inserted 

(fig. 4) between the radar and the command signals OC and ft 

L c *c 

when the system operates in the radar prime mode. It has also 
been shown (ref. 4) that if the radar operates in the prime mode, 
none of the radar components ^u R , v R or w R j can be absent 

for longer than a few msec or vehicle control will be lost. Thus, 
in the radar prime mode, the inertial navigator must be used as a 
backup to handle radar-dropout periods. In this case, equations 
(5) thru (10) are solved and the outputs are substituted for 

U R J v r 3 W R > and h a when the radar lock-status flags indicate 
the loss of a signal. 

Figure 9 shows the 6D0F simulation results for a system op- 
erating in the radar-aided inertial mode. Figure 10 shows the 
results for a system operating in the radar prime mode with in- 
ertial backup. Notice the difference in the effect of radar 
noise on the vehicle's thrust response. When the system is op- 
erating in the radar-aided inertial mode, it is about four times 
less sensitive to radar noise. 

Since the radar did not unlock on these runs, the radar prime 
mode did not use inertial navigator data except during the tip- 
up phase, when the attitude rates exceeded 10 deg/sec. This is 
a logic test in the simulation, since for high vehicle rates, 
the radar data lags significantly. Figures 9 and 10 show slight 
trajectory differences because, in each run, the navigator was 
initialized with errors. Even so because of radar data updating, 
these errors are reduced as time increases . 


Radar noise 
f il ter 



Figure 8.- Block Diagram of the Steering Channel Control-Loop Compensation Subsystem 
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Figure 9.- 6DOF Simulation Results; System Operating in Radar-Aided 
Inertial Mode 
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Figure 9.- Continued 
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Figure 9.- Concluded 
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Figure 10.- 6DOF Simulation Results; System Operating in Radar Prime Mode with 
Inertial Backup 
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Figure 10.- Continued 
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Figure 10.- Continued 
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A second alternative mechanization which we have considered 
involves using slant range of the radar instead of a broadbeam 
altimeter measurement for descent-contour guidance and for ini- 
tiating the terminal-descent phase. As far as descent-contour 
guidance is concerned, it makes little difference whether we use 
range or altitude, although a single range beam is somewhat more 
sensitive to surface— slope variations than the altimeter measure- 
ment. Figure 11 shows 6D0F simulation results for a radar-aided 
inertial system using range measurements. 

In this simulation, the Lander did not land successfully. 
Instead, because of the surface slope, the range had a large 
positive error. As a result, the Lander did not follow the de- 
scent contour during the early stages of the descent, and the 
thrust-to-weight ratio of the vehicle would have had to be in- 
creased for the vehicle to land successfully. 


C. Inertial Navigator Analysis 

Two methods of mechanizing the inertial navigator have been 
investigated. The first method, given by equations (5) thru 
(10), was based on a straight-forward solution of the vehicle 
translation equations of motion in the rotating-body-axis coordi- 
nate system. The errors in the outputs of the inertial system 
are then controlled by comparing them to the radar outputs and 
feeding back the differences through gains (Ku, Kv, etc) that 
are selected via simulation and a control stability analysis. 

The solution of equations (5) thru (10) requires three elements 
of a direction cosine matrix A, which relates the vehicle-body- 
axis coordinate system to the local vertical. This matrix solu- 
tion is mechanized in the vehicle's digital computer by integrat- 
ing the following equations: 


Ajli = 

(A21 

• r ) 

(A31 

’ 9 ) 

( 25 ) 

A 12 = 

(a 22 

* r) 

- (A32 

O 

( 26 ) 

A X 3 = 

(A23 

• r) 

" (A33 

9 ) 

( 27 ) 

A 21 = 

( A 3 l 

’ P) 

(An 

O 

( 28 ) 

A 22 = 

(A32 

• P) 

(A12 

0 

( 29 ) 

A23 = 

(A33 

• P) 

“ ( A i 3 

0 

( 30 ) 

A31 = 

(An 

* 9 ) 

" (A21 

p) 

( 31 ) 

A32 = 

(Al 2 

’ 9 ) 

“ (A22 

p) 

( 32 ) 

A33 = 
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' 9 ) 

- (A23 • 

p) 

( 33 ) 
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Figure 11.- Continued 
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Figure 11.- Continued 
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Figure 11.- Concluded 



The transformation equation is then 


x 


X_ 

b 


I 

Y b 

= [A] 

Y I 

z 



_ b_ 


I 


( 34 ) 


where X , Y^, and are the vehicle roll, pitch, and yaw 

axes, respectively, and X^ Yj, and Zj. are the local vertical 

coordinates ^Z^. points down^. 

Though only A l3 , A ss , and A^ are required in equations 

(5) thru (10), the digital computer will still be able to compute 
the entire A matrix. This is because equations (25 thru (33) 
must be solved during the deorbit phase to relate the body axes 
to the Sun-Canopus coordinate system of the orbiter. Thus, the 
A matrix computations must be reinitialized at some point so that 
they represent the transformation shown in equation (34). 

A good choice for this point is near the maximum-deceleration 
point during the entry. At this point, the angle-of -attack oscil- 
lations will be negligible, so that vehicle's roll axis will lie 
along the velocity vector. Figure 12 shows the range of values 
that the flight path angle can have at maximum deceleration. 7 M 

is the angle the velocity vector makes with respect to the local 
horizontal. If the vehicle is targeted for an entry angle, 7 E > 

of 16 deg, the 3a dispersion in this number can be -1 . 5 deg, as 
shown. The total dispersion of 7^ is then —3.7 deg (worst case). 

If the effect of the entry angle dispersion (-2.8 deg) is root-sum- 
squared with the effect of the atmosphere (-1.4 deg) the 3a dis- 
persion in 7 m i s ~3.1 deg. This error includes the effect of 

variations in the downrange distance due to variations in the atmos 
pheric density. It can be reduced in the future by three methods: 

1) Reducing the errors in that result from orbit 

determination errors; 

2) Refining the atmospheric model; 

3) Using inertial navigation to compensate for variations 
in the downrange angle. 
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Flight-path angle at maximum deceleration, deg 



Entry angle, 7 & , deg 


Figure 12.- Flight-Path Angle at Maximum Deceleration vs Entry Angle 



Further studies will define how much the 3.1~deg error can be 
reduced. Note that this error will occur only in the pitch plane. 
The yaw error will be much smaller. The method of updating the 
roll attitude is based on the fact that the vehicle* s yaw gyro 
rate will be zero if the pitch axis is horizontal. This is due 
to the normal gravity turn rate of the trajectory, which is as 
high as 1/2 deg/sec. 

In equations (5) thru (10), u, v, w, and H are initialized 
by being set equal to the first valid radar data that are received; 
then, the solution of the equations begins from that point. This is 
the simplest means of obtaining adequate performance from the iner- 
tial navigator. 

The second method considered for the mechanization of the in- 
ertial navigator was based on the use of a Kalman filter. This 
was done: 

1) To determine whether the more-complex, Kalman-f il ter 
approach was better from an operational standpoint; 

2) To determine whether the Kalman filter could be used 
to update the A matrix. 

A previous study on this subject (ref. 7) indicated that the 
Kalman-f ilter approach (see Appendix A) could produce safe land- 
ings with fewer than three Doppler radar beams operating. However, 
because that study was conducted for a nominal gravity-turn tra- 
jectory and did not consider the random effects of wind velocity, 
surface slope, and atmospheric density, the use of a nominal tra- 
jectory may be misleading. The approach used in this study was to 
mechanize a quasi-linear Kalman filter and simulate its operation 
in the M0D6MV program. This eliminates the assumption of a partic- 
ular nominal trajectory. 

The structure of this filter approach is shown in figure 13. 

An inspection of figure 13 shows that the only basic difference 
between using the Kalman filter and solving equations (5) thru 
(10) is that the Kalman filter generates a matrix of time-variable 
feedback gains, whereas a diagonal constant matrix (K , K , K , 

is used in equations (5) thru (10). We estimated that pro- 
graming the Kalman-filter approach would require about 500 addi- 
tional locations in the GCC . 

A further discussion of the modeling procedure is given in 
Appendix A. The results are briefly summarized on page 48. 
















1) The Kalman filter inherently initializes the inertial 
outputs with the first valid radar data. 

2) The Kalman filter appears to give more accurate esti- 
mates of u, v, w, and H during the first several 
sec of operation than the simple approach. 

3) Special care must be taken to keep the Kalman filter 
from ignoring the radar data as time progresses. If 
this is not done, there will be data saturation (un- 
detectable errors in the inertial estimates). 

4) More accurate inertial estimates do not produce a major 
payoff in terms of the descent performance, since the 
quality of the radar data improves with decreasing 
altitude anyway. 

5) The expansion of the Kalman filter from four states 
(u, v, w, and H) to seven (u, v, w, H, A 13 , 
A 23 , and A 33 ) was considered unnecessary in light 
of comment 4, particularly since doing so would have 
significantly increased the computer requirements. 

The errors of the inertial navigator can be estimated by refer- 
ring again to equations (5) thru (10). The error sources and the 
resulting acceleration errors are shown in table 1 . 


TABLE 1.- INERTIAL SYSTEM ACCELERATION ERRORS 


Error 

Source 

Error 

Equation 

Error 
Magnitude , 
ft/ sec 2 

Axis 

Accelerometer 
bias (400 pg) 

E = Bias 

0.005 

All 

Accelerometer 
Scale Factor 

(1/4%) 

E = S • A 

f X 

0.125 

Axial 

Gyro drift 
(1 deg/hr) 

E = u • Drift 

0.0024 

Lateral 

b 

Attitude error 
(10 deg) 

E = g • sin 10° 

2.1 

Lateral 

a u = 500 fps, A_^ = 50 ft/sec 2 ; p 
^Radar unlocked. 

, q, r, v, w, Ay, Az = 0. 



Table 1 shows that, in the absence of radar data, there will 
be significant errors in the lateral velocity channels (v and w) , 
and the side velocity will increase 2.1 ft/sec. These results were 
substantiated using the M0D6MV program. Thus, the inertial system 
should be considered only as a short-term memory device for extra- 
polating data during radar-dropout periods. 

When the radar is locked, the attitude error produces a lateral 
velocity bias. This can be derived from equations (5) thru (10) by 
considering A , p, and q to be zero and A 33 to equal the 

sine of the error angle e. 

w = g sin e + K w (w R - w) 

The transfer function of this equation is 

g sln 6 + ( K „) K) 

S + K 
w 

In the steady state (S = 0), we wish to obtain w = w„ which 

R 

will be true only if e = 0. When e = 10° and K =1, the 
steady-state error in w will be w 


E 

w 


g sin e 

K 

w 


2.1 fps. 


K > K „> K l 9 and K have been set equal to 1 in studies to 

U V w u 

date. Further studies should determine whether these gains can be' 
increased to decrease the effect of the attitude error. 


D. Radar Requirements 

Requirements must be established for the operating range, ac- 
curacy, probability of unlock, and probability of false lock be- 
fore selecting the final design of the radar. Analyses have been 
conducted to define these requirements. 


(35) 

(36) 

(37) 
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The point on the trajectory at which accuracy is most critical 
is where the engine cuts off just above the surface. The accuracy 
to which the cutoff conditions can be met affects not only the de- 
sign of the landing structure, but also the associated probability 
of a successful landing, as well as the amount of surface heating 
from the engines. Figure 14 shows the vertical impact velocity as 
a function of cutoff conditions. The shaded area is the range of 
conditions obtained using 30 cutoff conditions of: 

Altitude = 10 ft ± 5 ft; 

Velocity = 10 fps ± 3 fps. 

The currently specified values of the parameters important to 
landing success are shown in table 2. The expected performance is 
based on the performance of the modified LM radar. Other proposed 
designs would also be expected to meet these requirements. 


TABLE 2.- EFFECT OF RADAR ERRORS ON LANDING CONDITIONS 



Attitude , 
deg 

Attitude 
rat e , 
deg/sec 

Vertical 
velocity , 
fps 

Horizontal 
velocity , 
fps 

Altitude , 
ft 

Nominal engine 
cutoff 

0 

0 

10.0 

0 

10.0 

Nominal impact 

0 

0 

18.3 

0 

0 

Desired accuracy 
at cutoff (30) 


5.0 

3.0 

5.0 

5.0 

Expected perform- 
ance (30) 

H 

3.6 

0.7 

1.0 

2.5 

Desired accuracy 
at impact (30) 

6.0 

5.0 

5.0 

6.0 

-- 

Expected perform- 
ance (30) 

4.8 

3.6 

2.1 

1.0 



At higher altitudes, the accuracy is not as critical as it is 
near cutoff. The 30 high-altitude accuracy requirements are 
shown below: 


Altitude or 

slant range: ±4.5% or 5 ft (whichever is greater); 

Velocity: ±4,5% or 3 fps (whichever is greater). 
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Impact velocity, 





The maximum velocity that will be measured by the radar is 600 
fps. This velocity occurs during the parachute phase. Based on 
the above requirement, the 3 cr maximum velocity error is 27 fps, 
which is not significant when combined with random wind velocities 
on root -sum-squared (RSS) basis. Other M0D6MV simulations con- 
ducted during this study have shown that successful landings can 
be achieved with initial velocity errors up to 30 fps, and that 
the 3a altitude or range error due to radar tolerances of 4.5 % 
at the vernier ignition altitude of 4000 ft is 180 ft. Because 
successful landings have been simulated with initial errors of up 
to 150 ft, we feel that the specified accuracies are adequate for 
guidance purposes. 

The operating range of the radars depends on the overall tra- 
jectory of the vehicle and the design of the aerodynamic decelera- 
tor. The radars should be designed to overlap the staging alti- 
tudes to allow for some later adjustments in the trajectory. The 
LARA and the TDLR should operate between 20 000 and 10 ft for ve- 
locities ranging from 1000 fps to zero. The HARA should operate 
between 200 000 and 10 000 ft for velocities ranging from 20 000 
to 500 fps. 

The probability of obtaining valid radar data has been studied 
using Monte Carlo techniques. These will be described later in 
this report. An important question to be answered is: how long 

can we tolerate the loss of radar data? This question must be 
considered in terms of the alternative guidance modes described 
previously . 

Simulations conducted with the system operating in the radar 
prime mode have shown that the mission will be unsuccessful if 
one of the velocity components (one beam of a 3-beam radar or two 
beams of a 4-beam radar) unlocks during the tip-up maneuver (see 
ref. 4). Therefore, a suitable short-term backup mode is required. 
The best backup mode is to use is an inertial navigator. Without 
the navigator backing up the system, the radar would be required 
to have a 99.7% probability of maintaining the radar locks during 
the tip-up maneuver. 

The mission can also be a failure if the radar does not detect 
the altitude at which the vernier engines should ignite (4000 ft) 
during the parachute-descent phase. To preclude such a failure, 
the LM radar, when operating in the radar prime mode, should have 
a 99.7% probability of locking on the range beam and two of the 
velocity beams. 


If the system operates in the radar-aided inertial navigation 
mode, or in the radar prime mode with an inertial backup, these 
requirements can be relaxed. Since the inertial navigator is a 
short-term memory and data extrapolator , intermittent radar sam- 
ples can be tolerated. During the parachute descent, near vernier- 
ignition altitude, the vertical velocity of the vehicle will be 
between 100 and 300 fps. Then, if the 3a velocity accuracy is 
4.5 1 , the maximum velocity error will be 13.5 fps (table 3) and 
the drag acceleration will be about one Mars g. Table 3 shows 
how the altitude error and the velocity error will accumulate if 
radar data are lost for various times. The table shows that the 
axial velocity can be lost for up to 2 minutes with negligible 
effect. However, if the altitude data and lateral velocity data 
are lost for more than 10 sec, significant errors accumulate. 
Consequently, we recommend that, before the vernier engines ignite, 
at least one valid signal be received during the following intervals 

Altitude data: 10 sec or less; 

Lateral -velocity data: 10 sec or less; 

Axial-velocity data: 120 sec or less. 

These recommendations may be restated as follows: before the 

Lander descends to 4000 ft, valid radar data should be received 
more frequently than this; and, furthermore, the probability of 
receiving these data (i.e., the probability of locking the re- 
quired number of radar beams) should be at least 0.997. 

For the terminal -descent phase, the errors shown in table 3 
are essentially the same during the initial part of the trajectory. 
Only the accelerometer error is higher; this is because the meas- 
ured acceleration increases from about 1 g to about four g. 
However, if radar data are lost for less than 10 sec, the acceler- 
ometer error is still negligible. Since the terminal -descent phase 
lasts less than 40 sec in any case, it would be desirable to update 
the radar data more often than every 10 sec. 

Table 4 shows how the altitude error and the velocity error 
will accumulate as the result of radar losses during the constant- 
velocity descent phase. Using a 3a landing velocity of 6 fps, 
we see that a 3-sec loss of radar data in that channel is accept- 
able. However, using the specified requirements for altimeter and 
TDLR accuracy (see table 3), we see that the altimeter must provide 
data until the vernier engines shut down in order to hold the alti- 
tude accuracy to ±5 ft. If the axial-velocity data were accurate 
to 2 fps and the altitude to 2 ft, then the altimeter data could be 
lost for up to 2.5 sec. The axial velocity can be measured with 
this accuracy since it is the average calculated from at least two 
radar beams. The feasibility of limiting the range error to 2 ft 
depends on the mechanization that is chosen. 



TABLE 3.-3 a NAVIGATOR ERRORS DURING THE PARACHUTE-DESCENT PHASE 


Axial velocity 
a c J 
errors , fps 


Altitude errors, ft 


Lateral velocity 
errors , fps 


TDLR b 

Accel er- 
omet er 

RSS 

Altimeter* 

13.5 

0 

13.5 

180 

13.5 

.03 

13.5 

180 

13.5 

.15 

13.5 

180 

13.5 

.3 

13.5 

180 

13.5 

.9 

13.5 

180 

13.5 

1.8 

13.6 

180 

13.5 

3.6 

14.0 

180 


405.0 13.5 

810.0 54.0 
1620.0 216.0 


ometer RSS TDLR Attitude 

0 180 11 -2.1 

.015 180 11 0 

.37 191 11 8.4 

1.5 225 11 18.9 

13.5 443 11 60.9 

54.0 830 11 123.9 

216.0 1640 11 250.0 


Gravity and attitude errors negligible. 

^Vertical velocity = 300 fps, 3a radar accuracy = 4.5%. 

~ Axial acceleration = 12.3 ft/secf 3a accelerometer accuracy = 0.25%. 
3 

Altitude = 4000 ft, 3a altimeter or range accuracy = 4.5%. 
e Lateral velocity = 250 fps (worst case), 3a radar accuracy = 4.5%. 
^3a attitude error = 10°. 


TABLE 4.- 3a NAVIGATOR ERRORS DURING THE CONSTANT 
VELOCITY DESCENT PHASE 


Altitude errors, ft 


Radar unlock 
time, sec 


0 

1.0 

2.0 

3.0 

4.0 

5.0 


Altimeter 

TDLR 

RSS 

TDLR 

5.0 

0 

5.0 

3.0 

o 

in 

3.0 

5.8 

3.0 

5.0 

6.0 

7.8 

3.0 

5.0 

9.0 

10.3 

3.0 

5.0 

12.0 

13.0 

3.0 

5.0 

15.0 

15.8 

3.0 


Lateral v|locity 
errors, fps 


Attitude 


- 2.1 

0 

2.1 

4.2 

6.3 

8.4 



Altitude < 100 ft, velocity < 67 fps, axial -velocity error 
3 fps, gravity and accelerometer errors negligible. 


Attitude error 


navigator gain 





















From these studies, we determined that the radar data should 
be updated every three sec or less when the altitude of the Lander 
is less than 100 ft. When the altitude is between 100 and 4000 ft, 
the update interval, t ^ , can be calculated from equation (38). 

l ud - WO < H - 100 > + 3 <38) 

In there is a false lock, such as a main-lobe-to-cross -lobe 
lock or a side-lobe-to-side-lobe lock, the radar will put out 
false data that will be interpreted by the navigator as correct 
data. As documented in the next chapter, simulations have shown 
that a cross-lock condition which persists longer than three sec 
can cause loss of the vehicle. Thus, the probability of having a 
false lock for longer than 3 sec should be less than 0.003. 
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IV. MODIFIED LM RADAR STUDY 


Considerable effort was made to investigate the application 
of the Lunar Module (LM) radar and a modified version of it (the 
Mod-LM radar) for the Mars Lander mission. Among the tools used 
were the static Monte Carlo programs and the M0D6MV program. The 
results of these analyses are documented in references 4 thru 6; 
but, for convenience, this chapter contains a summary of the 
principal results of those simulations and a description of the 
results obtained using the time -correla ted Monte Carlo program. 

The geometry of the LM and Mod-LM radar beams is shown in 
figure 15. 
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Figure 15.- Beam Geometry for the LM and Modified LM Radar System 


The Mod-LM radar is a cw/FM-cw radar that uses three independent 
velocity trackers and one range tracker. Velocity compensation 
is required from Beams 1 and 3. The LM and Mod-LM radar para- 
meters are compared in table 5 . 



TABLE 5.- LM AND MODIFIED -LM RADAR PARAMETERS 


Parameter 

LM radar 

Mod -LM radar 

Transmitter power, dBW 

-13.0 

-13.0 

Losses, dB 

-1 .8 

-1 .8 

Receiver sensitivity 

a 

b 

Antenna gain, dB 

28.0 

28.0 

Maximum beam width, deg 

7 .5 

7.3 

Minimum beam width, deg 

3.7 

3.7 

Carrier frequency in velocity channel , GHz . 

10.51 

10.51 

Carrier frequency in range channel, GHz . . 

9.58 

9.58 

Sawtooth modulation frequency in range 

r 

c 

channel , Hz 

130.0 

130.0 

Upper high -mode search limit, KHz 

d 

e 12 . 0 

Lower high -mode search limit, KHz 

f 

e 0 

Upper low -mode search limit, KHz 

e 7.4 

e 6 . 0 

Lower low -mode search limit, KHz 

e -5 .24 

e 0 

Acquisition time, sec 

6 

3 

Mode -switching attitude, ft 

2500 

1250 

High -mode filter bandwidth, Hz 

2800 

600 

Low-mode filter bandwidth, Ha 

400 

400 


Velocity -dependent (3-dB signal -to -noise ratio used as criterion) 

^Velocity -dependent (3-dB s ignal -to -noise ratio used as criterion) 
More sensitive than LM radar. 

'Frequency-modulated . 

*52.0 KHz for Beams 1 and 2, 70.0 KHz for Beam 3. 

0 

For all four beams. 


-41.0 KHz for Beams 1 and 2, -20.0 KHz for Beam 3, 







The modifications made to the LM radar included adding a 
centered range beam, decreasing its search time, increasing its 
sensitivity, obtaining parallel data output on a per beam basis, 
and using a new range criterion for mode change. Note that before 
a production model of the Mod-LM radar can be heat-sterilized, 24 
of the 104 components and 51 of the 124 types of material that are 
not compatible must be changed. These changes were not included 
in the simulations of the Mod-LM radar. 

A. Summary of Previous Results 

The four techniques used to determine the application of the 
Mod-LM radar were: 

1) Tracking analyses (to determine the effect of changing 
the doppler frequency and the range); 

2) Error analyses (to qualify the noise models being 
used and to determine the effect of noise on the 
vehicle ) ; 

3) Monte Carlo analyses (to determine the probability of 
detecting various radar beams); and 

4) 6 DOF terminal -descent simulations (to determine noise, 
mode -switching , and system effects). 

The tracking rate capability of the Mod-LM radar is summarized 
in table 6. The only area of concern was the tendency of the range 
beam to unlock in the high mode during tip-up maneuvers. 


TABLE 6.- TRACKING RATE ANALYSIS RESULTS, MODIFIED 
LM RADAR 



Velocity beam 

Range beam 

Program 

High 

mode 

Low 

mode 

High 

mode 

Low 

mode 

Parachute 3 

OK 

Does not 
exis t 

OK 

Does not 
exist 

Tip-up 3 

OK 

Does not 
exis t 

Break 

i i b 

lock 

Does not 
exist 


Parachute and tip-up rates = 30 deg/sec. 


Break lock at 30 deg/sec occurs at an incidence angle 
of 60 deg (verified in 6D runs). 



To qualify the noise models being used, the M0D6MV program was 
’■frozen" at various points along the trajectory (indicated in the 
upper half of table V-3) and a large number of samples were taken 
of the tracker output and resulting body velocities. These data, 
in turn, were reduced statistically to probability density func- 
tions. The resulting roo t -mean-square (rms) levels were compared 
with data generated by Robert Harrington of Ryan Aeronautical Co, 
San Diego, California, at the same trajectory points. The results 
were the same (tables 7 and 8). The effect of radar noise on the 
vehicle may be seen in table 8. Note that when the system operates 
in the radar prime mode, the specified condition at engine cutoff 
is not met for attitude rate. 


TABLE 7.- ERROR ANALYSIS RESULTS 


Time 

30 radar error. 

f ps 

X-axis 

Y -axis 

Z-axis 

Vernier engine ignition 

1 .65 

7.2 

4.5 

Tipup 

1.1 

4.2 

3.0 

Contour point 

.9 

3.6 

1.8 

Vernier engine cutoff 

.2 

.8 

0.5 


TABLE 8.- 3c EFFECT OF ERRORS ON SYSTEM PARAMETERS AT 
VERNIER ENGINE CUTOFF 


Parame ter 

Operating mode 
of system 

Specified 
cutof f 
errors 

Radar 

Prime 

Mode 

Radar 

Inertial 

Mode 

Attitude rate, deg/sec 

12.6 

3 .6 

5.0 

Attitude , deg 

4.2 

1 .2 

1.0 

Horizontal landing velocity, fps 

1.0 

1 .0 

6.0 

Vertical landing velocity, fps . 

2.1 

2.1 

3.0 




Monte Carlo analyses were used to investigate the probabili- 
ties of beams unlocking, crosslobe lockups, and incorrect range 
for the Modified LM radar. The results of these analyses are 
shown in table 9 and in figures 16 and 17. The parametric con- 
ditions for the analyses were: altitude = 4000 ft; nominal pitch 

angle = 20 deg from the vertical, 3c ground slope = 34 deg, 3c 
wind velocity = 146 fps, and 3c attitude dispersion caused by 
wind gust = 22.4 deg. 

The probabilities of a single beam becoming unlocked or supply- 
ing incorrect range data were analyzed for the upper and lower 
surf ace -ref lec tivity models specified in reference 8 over the range 
of velocities caused by variations in the atmospheric models. 

TABLE 9.- CROSS-LOBE LOCK STATISTICS 


Conditions 

Total 

Three main 

Cross-lobe 

Cross -lobe 



conditions 
per beam 3 

V 

t 

°0 

cases 

lobe locks 

locks 

100 

Minimum 

1000 

606 

0 

0 

200 

Minimum 

1000 

933 

0 

0 

300 

Minimum 

1000 

978 

0 

0 

100 

Maximum 

1000 

988 

560 

4 

200 

Maximum 

1000 

1000 

1450 

0 

300 

Maximum 

1000 

1000 

2136 

0 


3 

A cross-lobe condition occurs when the cross-lobe power is 
above the detector threshold present in the beam. The 
tracker may or may not be locked to the main beam signal 
under this condition (see Chapter VI, Section A). 


Figure 16 shows that it is impractical to operate the existing 
or Modified-LM radar in a primary radar control mode if all three 
velocity beams are used on a continuous basis. However, it may 
be possible to operate the Modified-LM radar in the radar-aided, 
inertial -control mode, even though some risk of mission failure 
is implied by the high probability of beam unlock, particularly 
for the lower boundary of reflectivity characteristics and low 
ve locities . 

The probability of having range errors is so high, even for 
the Modified-LM radar, that it cannot be used to measure altitude 
for engine ignition and parachute release. This means that a 
separate widebeam altimeter or the range along each beam, is 
needed. The results of the cross-lobe lock study, which was made 
at various velocity conditions, are shown in table 9. 





e t = 2 



100 200 300 


Vertical velocity, fps 

Figure 17.- Monte Carlo Analysis Results: Range Error Probability, 

Modified LM Radar System 


Using a cross -coupling threshold of -28.9 dB (an average 
from five flight LM radars), the study shows that the conditional 
probability of cross -lock is 32 °L for = 100 fps, based on the 

maximum reflectivity curve. The conditional probability, P(c/u), 
is expressed in terms of the probability of a beam unlock, P(u), 
and the joint probability of a cross-lobe lock occurring simul- 
taneously with a beam unlock, P(c,u). A cross -lobe lock is 
catastrophic to the mission unless the condition can be detected 
and compensated for. Figures 18a thru 18d show that this should 
be done in less than three sec. This condition is hard to detect, 
because there are widely -vary ing dynamic changes on the parachute 
and capsule during the terminal -descent phase. 

Because of the major changes needed to use the LM radar for 
the Mars Lander mission and the uncertainties of operational 
risk, a new radar system should be developed. 

The machine plots shown in figures 18a thru 18d are examples 
of the type of evaluation that was done with the MOD6MV program. 
The Beam 2 radar signal was suppressed (forced) below its thres- 
hold to determine if the radar would cross -lock, and if so, how 
long the main beam radar signal would have to be below the thres- 
hold and the beam cross-locked to cause mission failure. By vary- 
ing the length of time the main beam signal was suppressed, it was 
shown that three sec is critical. The plots show that cross-lock- 
ing Beam 2 to a side to make its direction that of Beam 1, while 
suppressing the main beam signal for three sec caused the vehi- 
cle to crash with very high side velocities. The main beam sig- 
nal was suppressed eight sec into the terminal descent and 
allowed to return at 11 sec. The vehicle had been driven off 
so far, due to the error caused in velocity by the cross -lock, 
that the main beam spectrum was no longer in the tracker band- 
width, and the radar remained cross -locked until the vernier 
engines cut off. 

A large number of M0D6MV runs were made to determine the 
compatibility of the Mod-LM radar with respect to acquisition 
time, beam dropout, tracking rates, radar noise, mode switch, 
etc. Many of these results are documented in reference 4. 


Altitude, ft 








B. Time -Correlated Monte Carlo Analysis Results 

Table 10 summarizes the data obtained with the time -correlated 
Monte Carlo program described in Appendix B. The data are from 
four computer runs -- i.e., the tip-up and parachute programs run 
with the upper and the lower bounds of the reflectivity model de- 
fined in reference 8 (see fig. 23). Each block of data in table 
10 is numbered for reference in later discussions. 

Four types of probability data are presented: the data la- 

beled "length of time that three beams are locked" give the proba- 
bility that three beams will be locked or cross -locked for T 
sec somewhere on the trajectory (T is plotted from zero to the 
maximum time that occurred in the random trajectories computed); 
the data labeled "length of time ... before vernier engines ignite 
give the probability that beams will be locked or cross -locked 
before the vehicle descends to 4000 ft. (These data are plotted 
out only to the minimum trajectory time). 

Table 11 indicates the mean, maximum, and minimum durations of 
the trajectories sampled with the two programs. 

The data generated in figures 19 thru 21 are for 3-, 4-, and 
5-beam radars and are based on the MOD-LM radar parameters. The 
3 -beam system has a lamda pattern, the 4-beam system has a square 
pattern, and the 5 -beam system has four beams in a square pattern 
plus a centered beam. 

As one might expect, the probability that three beams will be 
locked increases with the number of beams being used. This can 
be seen by examining table 10. To determine what beam configura- 
tion is optimum, a tradeoff is required between the number of beams 
and the transmitter power. Figure 24 shows that the probability 
of locking 3 beams of a 3 -beam radar is less than the probability 
of locking three beams in a 4-beam system whose transmitter power 
is nine dB less. 

Because of the increased power required for a 3-beam system 
to have a lock-on probability equal to that of a 4-beam system, 
it is more reasonable to choose a 4-beam system. In addition, 
we feel that the G&C requirement can still be met using a 4-beam 
system . 



TABLE 10.- TIME -CORRELATED MONTE CARLO RESULTS 
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(a) Parachute Phase, Maximum Reflectivity Model 


Figure 19.- Time-Correlated Monte Carlo Analysis Results: Probability That 

Three Modified LM Radar Beams Will Be Locked for T sec 
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Probability that three beams will be locked for 



T, sec 


(b) Terminal Descent Phase, Maximum Reflectivity Model 


Figure 19 .- Continued 
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Probability that three beams will be locked for T sec 


99.99 



(c) Parachute Phase, Minimum Reflectivity Model 


Figure 19.- Continued 
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Probability that three beams vill be locked for T sec. 


99.99 



(d) Terminal Descent Phase, Minimum Reflectivity Model 
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Figure 19.- Concluded 










Probability that beam 1 is cross-locked for 
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(a) Parachute Phase, Maximum Reflectivity Model, Beam 1 


Figure 20.- Time-Correlated Monte Carlo Analysis Results: Probability That In- 

dividual Modified LM Radar Beams Will Be Cross-Locked for T sec 
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Probability that beam 2 is cross-locked for T sec 



(b) Parachute Phase, Maximum Reflectivity Model, Beam 2 
Figure 20.- Continued 
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Probability that beam 3 is cross-locked for T sec 



(c) Parachute Phase, Maximum Reflectivity Model, Beam 3 


Figure 20.- Continued 
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Probability that beam 4 is cross-locked for T sec 



Figure 20.- Concluded 
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Probability that three beams will be locked for T sec before the vernier engine ignities. 



Figure 21.- Time-Correlated Monte Carlo Analysis Results: Probability That Three 

Modified LM Radar Beams Will Be Locked for T sec Before the Vernier 
Engine Ignites 
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Probability that three beams will be locked for T sec before the vernier engine ignites 
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Figure 21.- Concluded 













Joint probability of a cross-lock per beam 



Figure 22.- Time-Correlated Monte Carlo Analysis Results: Probability That 

Individual Modified LM Radar Beams Will Be Cross-Locked for 
T sec Before The Vernier Engine Ignites 
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0 10 20 30 40 50 60 70 80 90 


Incidence , angle , deg 

Figure 23.- Terrain Reflectivity Model Defined in the 
Mars Engineering Specification 




Probability of locking three beams 



Transmitter power per beam, dB 


Probability of Locking All Beams from a 3-Beam Radar, % 


Figure 24.- Time Correlated Monte Carlo Analysis Results: Probability of 

Locking Three Modified LM Radar Beams vs Transmitter Power 
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TABLE 11.- TRAJECTORY TIMES 


Program 

Trajectory time, sec 

Mean 

Maximum 

Minimum 

Terminal descent 

33 

45 

27 

Parachute 

60 

97 

37 


For the parachute descent and the maximum reflectivity model, 
all beams are locked or cross-locked at an altitude of 4000 ft 
(see Blocks 3 and 4 of table 10). However, the G&C requirement 
of one valid 3 -beam update every 10 sec (see Section III-D) can- 
not be met with either the 3-, 4-, or 3-beam radar, because of 
the high number of cross -locks. At first glance, it would seem 
that the 5-beam system for which there is a 1007> probability 
that three beams will be locked for 37 sec before reaching the 
altitude at which the descent engine ignites, would be adequate; 
however, there is no practical means for detecting which of the 
beams is providing good data (i.e., it is not feasible to detect 
which beams are cross -locked) . 

From Block 11 of table 10, we see that the probability of 
locking three beams before the Lander descends to 4000 ft is 
very low; only the 5-beam radar system gives a reasonable proba- 
bility. Note, however, that there are no cross -locks (Block 12). 

With a 4-beam radar we are 2.7% more likely to lock three 
beams at 4000 ft, and 217o more likely to lock three beams 37 sec 
before the Lander reaches 4000 ft, than if we use a 3 -beam radar 
sys tern. 

If we use the parachute descent program and the minimum re- 
flectivity model (Block 1 of table 10), we see that there is a 
99 . 27o, 99.7%, and 1007o probability of having three beams locked 
for at least 30 sec for the 3-, 4-, and 5-beam radars respectively. 
For the minimum reflectivity model (Block 9), the probabilities 
of locking three beams for the same length of time (30 sec) are 
70.77* for the 3 -beam radar, 897o for the 4-beam radar and 1007, 
for the 5 -beam radar. Only with the 5-beam radar are three beams 
locked for a satisfactory length of time. 



During the terminal -descent phase, all radars will have, three 
beams locked with 100% probability for the 27 sec before the en- 
gine cuts off, for both the maximum and minimum reflectivity 
models. The probabilities for the total time that three beams 
will be locked for each type of radar system are indicated in 
Blocks 5 and 13 . 

Actual cross -locks occurred during the parachute descent, 
using the maximum reflectivity model (Blocks 2 and 4). Most of 
these occurred for only a short time, such as when a wind gust 
tipped the vehicle and caused a beam to unlock and then cross - 
lock. However, a few times (two trajec tor ies /beam) , the beams 
cross -locked during the initial search and remained cross -locked 
until the Lander descended to 4000 ft. As mentioned before, this 
situation is not tolorable. 

In summary then: 

1) With the reflectivity models of figure 23, the MOD-LM 
radar gave satisfactory probabilities both for having 
three beams locked and for having no cross -locks dur- 
ing the terminal descent . 

2) During the parachute descent, using the maximum re- 
flectivity curves, the 4- and 5-beam radars gave 
satisfactory probabilities of locking three beams; 
however, for both systems, there were an intolerable 
number of cross -locks . 

3) The results obtained using the minimum reflectivity 
curve indicate that none of the radars performed 
satisfactorily. Accordingly, we made another time- 
correlated Monte Carlo run using a reflectivity model 
supplied by Dr. Richard F. Broderick (see figure 25). 
The upper -bound curve shown in this figure produces 
the greatest number of unlocks and was used in this 
new run. The radar that was simulated was a cw radar 
similar to the MOD-LM radar. The results of this 
analysis (figures 26 and 27) show that for the 4-beam 
radar, the probability of locking three beams satisfied 
the G&C requirement, and that, for the 5 -beam radar, 
the probability of having three beams locked for the 

37 sec before the Lander descended to 4000 ft was 100%. 







Probability that three beams will be locked for T sec 



f ill 1 Ml I ] 1 i M 11 1 m i II 1.1 1 1 Ml i i.Liu m II m Mill ill 11 11 II i m j j 11 II 11 II m II m 1 1 

0 10 20 30 40 30 60 70 80 90 


T, sec 

Parachute phase, maximum reflectivity model from figure 25 

Figure 26.- Second Time -Correlated Monte Carlo Analysis Results: Probability That Three 

J Radar Beams Will be Locked for T sec 
o 
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Probability that three beams will be locked for T sec before the vernier engine ignites 



0 5 10 15 20 25 30 35 40 45 

T, sec 


Parachute phase, maximum reflectivity model from figure 25 

Figure 27.- Second Time -Correlated Monte Carlo Analysis Results: 

Probability that Three J Radar Beams will be Locked 

o 

for X sec Before the Vernier Engine Ignites 
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Most cross-locks existed for a very short time; how- 
ever, when they occurred during the initial acquisi- 
tion, they were apparently maintained throughout the 
terminal descent. Such cross-lock conditions are in- 
tolerable . 

The 4-beam radar that had a square beam pattern will 
meet the G&C requirements and appears to be the best 
choice in terms of the tradeoff of transmitter power 
and the probability of having three beams locked. 


V. BESSEL SIDEBAND RADAR STUDY 


A technique to achieve isolation from receiver mixer noise (a 
consideration for the velocity sensor) is frequency modulation of 
the carrier. The following is an evaluation of a Bessel sideband 
radar in which data is processed at the J x sideband. 

The radar is a modified AN/APN-187 Doppler radar system. The 
general characteristics of this radar are: a 4-beam sequential 

lobing system; a time-shared sine-cosine tracker; Bessel-.^ side- 
band data processing; a program bandwidth for the tracking filter; 
and sequential frequency modulation on two frequencies . 


A, Radar Math Modeling 

A block diagram of the J! radar circuitry that was simulated 
in the M0D6MV is shown in figure 28. The following is a descrip- 
tion of the detailed parameters and calculations contained in each 
block . 

-L: Beam switching (Block A) .- The beam-switching rate of the 

radar is 24 Hz. This rate is maintained in the simulation by a 
counter that switches beam positions every 1/96 sec. The radar as 
modeled uses the MDD-LM velocity-beam geometry shown in figure 15 

2. Modulation frequency switching logic (Block B) .- The trans- 
mitter sequentially switches between 148 and 161 .5 KHz at a rate 
of 12 Hz. A counter is used to switch these modulation frequen- 
cies every 1/12 sec. 

3 . Bessel coefficient (Block C) .- The Bessel function gener- 
ator used is a subroutine called BESJ taken from reference 9. The 
Doppler processing is done at the frequency of the sideband. 

The argument for the Bessel function is 


2(1) sin (2 jtR/A) 

where the deviation index I is 1, the range in the beam is R, 
and the wavelength of the carrier A is 0.074 ft. 

Geometry (Block D) .- The selected ground slope is read into 

the program and is used to define a unit vector perpendicular 

to the terrain. Then, the incidence angle of each beam is calcu- 
lated in the following manner: 


(39) 
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Compute tracker 
bandwidth 














(40) 


J^Ij] ~ [ ^ ]T (^B j) ’ where J = 1> 2, 3 

C ° S 6 j = M Ij ‘ (41) 

where is the pointing direction of the jth beam in inertial 

coordinates, [A] is the body- to inertial-coordinate-system 
transformation matrix, is the jth beam to body-coordinate- 

system transformation matrix, ^ is the jth beam unit vector 

expressed in body coordinates, and 6 is the incidence angle as- 
sociated with the center beam propagation vector. 


The range along each beam is obtained from equations (42) and 
(43) . 



R = D/cos 6; 

(42) 


D = X • p 

n 

(43) 

where X is the position vector of the 
perpendicular distance from the vehicle 

vehicle and D is the 
to the terrain . 

The Doppler frequencies of the beam are calculated by: (1) 

transforming the inertial velocities of the vehicle into body-axis 
velocities; (2) transforming the body-axis velocities into beam 
velocity components; and (3) converting beam velocities into Dop- 
pler frequencies . 


V B . [A] Vj 

(44) 


V B1 - 4*1 • \ 

(45) 


V B2 = |J B2 ’ 4 

(46) 


V B3 = ^B3 ’ 4 

(47) 


where [A] is_the inertial- to body-coordinate-system transforma- 
tion matrix, V^. is the vector that expresses the inertial veloc- 
ity of the vehicle, V g is the velocity of the vehicle expressed 

in body coordinates, and V_ 1# V , and V_ 0 are the velocities 

151 Jdz dj 

of Beams 1, 2, and 3, respectively. 
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The Doppler frequencies are: 


f , . = 2V . / A , . < 

dj Bj/ dj 

where f is the Doppler frequency in a beam and A , . is the 
d j a J 

wavelength of that beam. 

5. Reflectivity model (Block E) The reflectivity model used 
for this M0D6MV radar simulation is shown in figure 23 . The two 
bounds shown in the figure include all the reflectivity models that 
have been proposed for the lunar surface. A curve-fit was obtained 
for each bound (the respective equations are keyed on input) . 

6. Signal bandwidth (Block F) .- The spectral bandwidth of a 
Doppler signal is calculated as follows: 


/\f = -^21 v 

^ d A B 


sin y 


where Ar is the average 3-dB beamwidth of the antenna (0.069 ra- 
dians) , V is the maximum velocity of the vehicle relative to 
B 

the reflecting surface, and y is the angle between the center- 
line of the beam and the velocity vector. 


Thus, the spectrum bandwidth is proportional to the component 
of velocity perpendicular to the beam-pointing direction, as modi- 
fied by the 3-dB beamwidth of the antenna. The spectrum is assumed 
to be Gaussian. 


7, Tracker bandwidth (Block G) .- The tracker of the Jj. radar 
contains filters whose bandwidths depend on the mode in which the 
radar is operating. The radar can operate in three different modes: 

Mode 1 -- Initial search mode in which the filters have a band- 
width equivalent to the maximum expected Doppler fre- 
quency (16 000 Hz) 

Mode 2 -- Final phase of search during which the filters have 
a bandwidth 1/8 of that in Mode 1 (2000 Hz) 


Mode 3 -- Normal track mode during which the filters have a 
bandwidth which is proportional to the velocity 
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When the radar operates in Mode 3, we have assumed that 

the velocity vector is always centered in the beam pattern. The 
filter bandwidth is then programed as the product of a coefficient 
times the total velocity of the vehicle. The coefficient is 

C = 2Ay sin y (50) 

where Ay is the beamwidth of the antenna (4 deg) and y is the 
beam splay angle (20 deg) . 

8. Signal- to -noise ratio (Block H) . - Two signal - to-noise (S/N) 
ratios are computed: one for the tracker filters and the other 

for the threshold detector. The latter has a fixed bandwidth of 
100 Hz. The S/N ratio may be defined in a peak-power, spectral- 
density sense. From figure 29, the S/N ratio according to this 
definition would be SI/NO. This definition was maintained in the 
simulation , 



Frequency 


Figure 29.- Power Spectral Density Used in the J ± Radar Simulation 
The received power may be calculated from equation (51) . 

P A 2 G W a L E 

pn — 2 

(4it) 2 R 2 cos G 

where PR is the power received, P is the transmitter power 

(167 mW) , A is the carrier- frequency wavelength (13.325 GHz), 

G is the antenna gain (30 dB) , W is the antenna weighting fac- 
tor (3.5 dB) , a is the relative cross-section of the radar, 
o 

L is the total microwave loss (5.1 dB) , and E is the Bessel 
coefficient. The noise spectral density (dB/Hz) is calculated 
from equation (52) . 
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( 52 ) 


where K is Boltzmann’s constant, T = 290 °K, and N is the 
receiver noise (9 dB) . 

Next, the received power is numerically integrated for the 
tracker filters and for the detector. As demonstrated in figure 
30, the signal spectrum and filter responses ideally have rectan- 
gular shapes. Power integration is required when the tracking 
rates induce a lag between the output of the voltage control oscil- 
lator (VCO) and the incoming Doppler signal . 


abw— 

• ) 


Center of 

Doppler 

spectrum 



Center frequency of 
tracker filter 


Figure 30.- Idealization of Received Power 


When this situation is detected in the simulation, the received 
power is modified as shown in equations (53) and (54) . 


PrT = PR + 10 log 



/ ABW D 

P RD * PR + 10 108 


where P is the signal power in the tracking filter bandwidth, 
AB W t is the proportional bandwidth of the signal in the tracking 
filter (see fig. 30), Af rf is the signal bandwidth, P RD is the 
signal power in the detector, and ABW Q is the proportional band- 
width of the signal in the detector (see fig. 41) . 

To obtain an approximate estimate of the peak power level of 
the signal spectrum in the tracker and in the detector, the powers 


(53) 

(54) 
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P and V are divided by ABW to obtain dB/Hz. Then, the 
K1 KU t 

S/N ratio can be obtained from equation (52) . 


9. Fluctuation noise (Block I) .- The statistical character- 
istics of the tracker output are obtained by adding the appropriate 
noise to the Doppler frequency f . Because the signal is assumed 

to have a Gaussian distribution, the rms level of the frequency 
jitter o and the correlation time of the noise T are functions 
of the Doppler spectral bandwidth. 



(55) 


JL 

irAf J 

d 


(56) 


10. Bias (Block J) Although we could have added biases to 
the signal, we did not; the outcome of the simulation was predicta- 
ble and nothing would have been gained by including them. 

11. Tracker (Block K) As shown in equation (57) , the closed- 

loop transfer function of the tracker was represented by a 

first-order lag. 

1/(1 + ts) (57) 


where T = 0 .3 sec . 

12, Threshold detector (Block L) The threshold detector was 
modeled to allow a different criterion (input variable) for each 
mode; 3 dB is the criterion for normal tracking. 

13. Mode and search logic (Block M) If a beam becomes un- 
locked on the proposed system, all beams go into Mode 1. During 
this search, the bandwidth of the tracker filter is held wide open 
for one sec. If, at the end of one sec, all beams are above the 
threshold, then the radar goes to Mode 2; if not, the radar is 
held in Mode 1 for another sec. If the radar is operating in Mode 
2 and at any time during a 1-sec interval all the beams are above 
the threshold value, the tracker immediately switches to Mode 3; 
but, if at the end of one second this has not occurred, then the 
tracker reverts to Mode 1 (see table 12) . 

The only signal driving the VCO is the error signal of the 
discriminator (i.e., when the radar is first turned on, the VCO 
will be randomly positioned, and is not driven in a search mode). 
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TABLE 12.- TRACKER MODES 


Mode 

F unc tion 

Tracker filter 
bandwidth, Hz 

Time, sec 

i 

Initial search 

16 000 

i 

2 

Search 

2 000 

<i 

3 

Normal track 

Variable 



The flags that indicate "data good" are also provided by the logic. 

B. MOD6MV Results 

The J x radar, as it was modeled here, was not found to be ac- 
ceptable. The simulation showed that: 

1) The bandwidth of the detector was too narrow (100 Hz) ; 
break-lock conditions frequently existed due to noise 
spikes and tracker lags . 

2) The variable-bandwidth filter in the tracker does not 
work well; the coefficient cannot be determined unless 
we assume an angular distance between the velocity 
vector and the centerline of the beam exists. The 
angle assumed must be equal to the beam splay angle; 
but, even so, this causes the filter bandwidth to be 
too narrow whenever the roll axis of the vehicle is 
not aligned with the velocity vector. 

3) The detector has no memory (i.e., momentary fading of 
the signal causes the radar to break-lock immediately) . 

4) The tracker is not allowed to break-lock on one beam. 
Consequently, if the signal to a beam is lost, all 
channels must revert to Mode 1 . 

Figures 31a thru 31c are shown in order to clarify the above 
conclusions. These curves were obtained from computer runs of the 
M0D6MV. The shaded areas in the figures represent the times at 
which the beams were locked. This run imposed very moderate con- 
ditions -- no ground slope, and a low velocity (initially 150 fps) 
in a vertical descent from 4000 ft -- and there should have been 
no tendency for the radar to break-lock. Nonetheless, numerous 
break-locks did occur, because the noise and acceleration momen- 
tarily drove the signal out of the tracker and detector bandwidths . 
Since the tracker had no memory, break-locks on all beams occurred 
immediately . 
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U-channel velocity, fps 



Figure 31.- M0D6MV Computer Program Results: Effect of Beam Locking 

on Velocity, Radar System 


I 


99 








Probability of data loss. 


The sta 
used to det 
the J x ra 
considered 
in figure 3 


H « 4000 ft 
attitude » ! 
upper reflei 
model 


Vertical velocity, fps 


Figure 32.- Static Monte Carlo Analysis Results: Probability 

of Data Loss for the Radar System 

The reflectivity model used is the same as shown in figure 23 , 
The figure indicates that there is a high probability of losing 
data. Some of the reasons for this are: 

1) A range hole occurs at 6400 ft (see fig. 33); 

2) As a result of sequential lobing, the efficiency factor 






Power coefficient 


3) 


Each time a channel unlocks, all beams are declared 
unlocked ; 

4) Whenever a break-lock occurs with one-frequency modu- 
lation, an additional efficiency factor of -3 dB is 
used . 

In contrast to these results, using the minimum reflectivity model 
resulted in a loss of data 1007, of the time. 

Figure 33 shows that, during the descent, the J x radar has 
a range hole at about 3000 ft. This was probably responsible for 
some of the unlocks shown in figure 31 . 



Figure 33.- Bessel Power Coefficients for the J x Radar System 

The preceding Monte Carlo and M0D6MV results are presented in 
response to the contract task to model and analyze a Bessel side- 
band radar. The particular design used in this evaluation turned 
out to be poorly suited to the mission. However, certain features 
and approaches which were used are attractive, and further evalua- 
tion should be made with better-chosen parameters. Unfortunately, 
time did not permit an iterative design analysis during this study. 
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VI. ALTERNATIVE RADAR MECHANIZATION STUDY 


A. Radar Modeling and Approach 

The static Monte Carlo program described in Appendix C was 
used as a tool for obtaining tradeoff data with regard to modula- 
tion techniques and cross -coupling . Three modulation techniques 
can be considered for the TDLR. These are: (1) continuous wave 

|cw-J Q j modulation; (2) Bessel sideband (Ji) modulation; and (3) 

interrupted continuous wave (ICW) modulation. To compare the 
three modulation techniques, three static Monte Carlo radar sub- 
routines were created, each representing one type of radar. Each 
radar was modeled in such a way that the parameters which were 
not dependent on modulation technique were equivalent. 

The following were equal for each radar: 

Antenna gain = 25 dB 
Transmitter power = -16 dBw 
Carrier frequency = 13.33 GHz 
Total losses = -8 dB 
Noise figure = 11 dB 
Tracker bandwidth = 600 Hz 

In addition, an S/N criterion of 3 dB was used for determining 
a lock condition (for each of the three systems) . 

The bases for tradeoff were: the probability of locking three 

beams; the probability of cross-locking; the probability of con- 
ditional cross-locking; the probability of side lobe locking; and 
the probability of conditional side lobe locking. 

Cross-locking refers to locking a channel to the side lobe 
signal of a receiver that coincides with the main beam of a 
transmitter. Conditional cross-locking occurs when a side lobe 
signal of a receiver is above a threshold value and the channel 
is locked to the main beam signal (again the side lobe of a re- 
ceiver coincides with the main lobe of a transmitter). Side 
lobe locking occurs when a channel is locked to a received side 
lobe signal and the side lobe of a receiver coincides with the 
side lobe of a transmitter. Conditional side lobe locking occurs 
when received energy at the side lobe is above the threshold and 
the channel is locked to the main beam signal. 



The tradeoff criteria will be to determine, for each modula- 
tion technique, the side lobe levels that will prevent side lobe 
locking at the vernier-engine-ignition altitude and the side lobe 
level (or isolation due to frequency diversity) that will be re- 
quired to prevent cross -locking at the vernier -engine- ignit ion 
altitude . 

The J q radar was used as a baseline for establishing the 

S/N criteria on the receiver noise figure and the tracker band- 
width. Note the receiver sensitivity curve shown in figure 34. 

A curve-fit for each of the curves of figure 34 was obtained and 

used as a detector model for the J radar. Since these curves 

o 

are based on a 3-dB S/N threshold, that threshold was used for 
the other two radars as well. The wide-band mode of the J 

o 

system specifies a tracker filter bandwidth of 600 Hz, Since 
this mode includes the major portion of the expected range at 
Doppler frequencies, this bandwidth (BW) was used consistently 
to calculate the tracker noise level in the other two systems. 

Note that if the wide-band -mode receiver sensitivity curve is 
extrapolated to higher frequencies, it would probably level off 
at approximately -162 dB. Figure 34 also indicates that the noise 
level is about -165 dB. 

From this, then, we can use equation (58) to calculate an 
equivalent noise figure at high frequencies for the J q system. 

N = KTBN (58) 

Eq 

where K is Boltzman's constant, T = 290°, B is the bandwidth 
of the tracker, and N is the receiver noise. The noise is 11 
dB. This value was used to calculate the noise power for the 
ICW and J ± systems. 

Each modulation technique forces a unique degradation to the 
required signal power for normal acquisition and tracking. Due 
to mixer noise in the J system, the signal power must be ex- 
pressed as a function of velocity (see fig. 34). The ICW system 
suffers an efficiency loss of approximately 4 dB, due to the duty 
cycle and receiver blanking. The system, besides operating 

on reduced sideband power, is range -sens it ive (see fig. 33), so 
the parameters used for the J j_ modulation were a deviation index 
of one and sequential switching of two modulation frequencies 
(148 and 161.5 KHz); in addition, if the power received using one 
of the modulation frequencies resulted in a signal below the de- 9 
tector threshold, then the power received using the remaining modu- 3 
lation frequency was degraded by 3 dB, due to the 50% duty cycle 
that would result. 
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Figure 34.- Acquisition Sensitivity of the J Radar Receiver 




If a sequential lobing radar had been used it would have had 
an efficiency loss of 6 dB; because of this we^assumed that such 
a radar would transmit 6 dB more power. Hence", the results that 
will be shown here for the J x system are valid for either a 
time-dependent or simultaneous lobing system. 

The reflectivity models used for the static Monte Carlo runs 
described here are shown in figure 25 . These models were furnished 
by Dr. Richard F. Broderick, and are felt to be representative of 
the radar cross-section for the 4000-ft altitude at which these 
computer runs were made. 

The Monte Carlo program provides the following types of data. 
(1) the probabilities of having one, two, three, and four beams 
unlocked (equivalent to the probabilities of having four, three, 
two, and one beam locked); (2) the conditional probability of 
having a cross-lock (i.e., the probability that the received power 
from a side lobe is above the threshold); and, (3) the joint prob- 
ability of having a cross -lock (the probability that a main beam 
is unlocked at the same time that a cross-lock occurs). 

In determining the probabilities of having locks and condi- 
tional locks, data were obtained for a 4-beam radar whose beams 
were separated by 90 deg and whose splay angles were 20 deg. The 
side lobe-side lobe data were obtained for a main beam splayed at 
20 deg with three side lobes in a plane defined by the centerline 
of the main beam and the roll axis of the vehicle, as shown in 
figure 35. 
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Figure 35,- Side Lobe Configuration Used in the Static Monte 
Carlo Tradeoff Analysis 
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The data obtained from this configuration are applicable whenever 
a side lobe of the main beam lies within a volume formed by ro- 
tating the center line of side lobe 1 about the centerline of the 
vehicle. This is so because the Monte Carlo program selects 
random values for the roll of the vehicle, the ground slope 
azimuth, and the wind azimuth. For example, the statistics for 
a side lobe of the main beam positioned at 4 in figure 35 will 
be the same as those for side lobe 2. The results of this in- 
vestigation will be described in the following section. 


B. Comparative Monte Carlo Results 

The results of the static Monte Carlo computer runs are shown 
in figures 36 and 37 . The radar parameters and reflectivity 
models that were used in these runs have been described in the 
previous section. The random variables are shown in table 13. 

All runs were made at an altitude of 4000 ft. 


TABLE 13.- RANDOM VARIABLES USED IN THE 
STATIC MONTE CARLO TRADEOFF ANALYSIS 


Parameter 

Distribution 

3a Value 

Horizontal wind 
Ground slope 

Attitude due to wind gusts 
Vehicle roll 
Ground slope azimuth 
^Wind azimuth 

Gaussian 
Gaussian 
Gaussian 
Uniform + 180° 
Uniform + 180° 
Uniform + 180° 

146 fps 
34 deg 
22.4 deg 


To determine the probabilities of cross -coupling for the 
various side lobe-to-main lobe combinations, the following param- 
eters were used: upper and lower reflectivity models (see fig. 

25); terminal velocity and nominal attitude paired as 100 fps and 
0 deg, 200 fps and 10 deg, and 300 fps and 20 deg; and the tip- 
up and parachute Monte Carlo programs. In the past, the terminal 
velocity and the nominal attitude had been paired as 100 fps and 
20 deg, 200 fps and 20 deg, and 300 fps and 20 deg; but in using 
the time-correlated Monte Carlo program we noted that the first 
two paired conditions do not exist. 

The results of this run are shown in figures 36a thru 36 i. 

The side lobe-side lobe data are shown in figures 37a thru 37f, 
and were obtained using the upper reflectivity model, a terminal 
velocity of 300 fps, and a nominal attitude of 20 deg for both 
the tip-up and the parachute programs. 




One-way side lobe level, dB 

(a) Parachute Program, Lower Reflectivity Model, Terminal Velocity = 300 fps, 
Nominal Pitch Attitude = 20 deg 

Figure 36.- Static Monte Carlo Tradeoff Analysis Results: Conditional Probability 

of Cross -Lock at Vernier Engine Ignition 


Conditional probability of cross-lock, 



(b) Parachute Program, Lower Reflectivity Model, Terminal Velocity - 200 fps. 
Nominal Pitch Attitude = 10 deg 

Figure 36.- Continued 
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One-way side lobe level, dB 


(c) Parachute Program, Lower Reflectivity Model, Terminal Velocity = 100 fps 
Nominal Pitch Attitude = 0 deg 


Figure 36.- Continued 



Conditional pi 





One-way side lobe level, dB 


(e) Tip-Up Program, Lower Reflectivity Model, Terminal Velocity = 200 fps 
Nominal Pitch Attitude - 10 deg 


Figure 36.- Continued 



Conditional probability of cross-lock. 


100 



One-way side lobe level, dB 

(f) Tip-Up Program, Lower Reflectivity Model, Terminal Velocity - 100 fps, 
Nominal Pitch Attitude = 0 deg 


Figure 36.- Continued 
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Conditional probability of cross-lock, 



One-way side lobe level, dB 


(h) Parachute Program Upper Reflectivity Model 
Nominal Pitch Attitude = 10 deg 


Terminal Velocity = 200 fps. 


Figure 36.- Continued 
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One-way side lobe level, dB 


( i) Parachute Program Upper Reflectivity Model, Terminal Velocity - 100 fps 
Nominal Pitch Attitude = 0 deg 


Figure 36.- Continued 




( j) Tip-Up Program, Upper Reflectivity Model, Terminal Velocity = 300 fps. 
Nominal Pitch Attitude = 20 deg 


Figure 36.- Continued 
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One-way side lobe level, dB 

(k) Tip-Up Program, Upper Reflectivity Model, Terminal Velocity = 200 fps 
Nominal Pitch Attitude - 10 deg 


Figure 36.- Continued 




Conditional probability of cross-lock, 


100 



One-way side lobe level, dB 

(1) Tip-Up Program, Upper Reflectivity Model, Terminal Velocity = 100 fps 
Nominal Pitch Attitude = 0 deg * 

Figure 36.- Concluded 
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Conditional probability of side lobe lock, 



Two-way side lobe level, dB 

(a) Parachute Program, Side Lobe 20 deg from Main Beam 

Figure 37. ~ Static Monte Carlo Analysis Results: Conditional Probability 

of Side-Lobe Lock at Vernier Engine Ignition, Calculated Using 
Upper Reflectivity Model with Terminal Velocity ** 300 fps and 
Nominal Pitch Attitude at 20 deg 







Conditional probability of cross-lock, 
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Table 14 indicates the results of the first part of the modu- 
lation tradeoff (i.e., the probability of having three and four 
beams locked for a 4-beam radar) . Note that the different modula- 
tion techniques produce nearly the same results. It appears that 
the most favorable radar is the ICW, the next most favorable is 

Ji radar, and the least favorable radar is the J . It should 
x o 

be pointed out that some data (not shown) were taken using a trans- 
mitter power of -10 dB for each radar, and there was even less 
perceptible difference between the modulation techniques at that 
power level. The beam unlocks are due to large incidence angles 
(ICW) , large incidence angles and low velocities (j , and 

large incidence angles and range holes (Ji). 

TABLE 14.- PROBABILITY OF FOUR BEAMS LOCKED/ PROBABILITY OF THREE 
BEAMS LOCKED FOR 4 -BEAM RADARS IN PERCENT 


Program 

Model 

V T , fps 

0 T > deg 

Modulation technique 

Ji 

J 

0 

ICW 

Tip-up 

Upper 

100 

0 

99.8/100 

99.9/100 

99.9/100 


reflectivity 

200 

10 

100/100 

100/100 

100/100 



300 

20 

100/100 

100/100 

100/100 

Tip-up 

Lower 

100 

0 

98.8/100 

98.9/100 

99.7/100 


reflectivity 

200 

10 

100/100 

100/100 

100/100 



300 

20 

100/100 

100/100 

100/100 

Parachute 

Upper 

100 

0 

99.8/100 

99.3/100 

99.9/100 


reflectivity 

200 

10 

100/100 

100/100 

100/100 



300 

20 

100/100 

100/100 

100/100 

Parachute 

Lower 

100 

0 

99.8/100 

98.9/100 

99.9/100 


reflectivity 

200 

10 

100/100 

99.9/100 

100/100 



300 

20 

100/100 

100/100 

100/100 


Table 15 indicates that the J modulation technique is more 

o 

susceptible to cross-locks. This implies that conditions under 

which the J system unlocks and cross-locks (side lobe-to-main 
o 

lobe) are zero Doppler conditions, to which the other radars are 
not susceptible. 
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TABLE 15.- NUMBER OF CROSS LOCKS AS A FUNCTION OF SIDE LOBE LEVEL 3 


One-way side lobe level, dB 

Modulation technique 

Ji 

J 

O 

ICW 

20 

i 

10 

1 

25 

i 

9 

1 

30 

i 

6 

1 

35 

i 

3 

1 

37.5 

0 

2 

0 


Terminal velocity = 100 fps, nominal attitude = 0 deg, upper 
reflectivity model, parachute program. 


On the other hand, figures 36a thru 362 show that, in general, 
the J system is less susceptible to conditions of cross-lock 

than either the or ICW. The data indicate that if the one- 

way side lobe level is approximately 43 dB, there is virtually no 
possibility of having a cross-lock on any of the radars (for the 
parameters used in this analysis). Although it is difficult to 
achieve a one-way side lobe level of 45 dB in the direction of 
the main beams of the transmitter for a 4-beam system (the mean 
for 5-beam LM radars is 28.9 dB) , each transmitter could be set 
at a different carrier frequency, so that the isolation would be 
greater than 80 dB. 

Figures 37a thru 37f show the probabilities of having con- 
ditional side lobe locks. Although the difference in suscepti- 
bility here is practically insignificant, the Jj_ system is the 
least susceptible and the ICW is the most susceptible. Note 
that side lobes with small splay angles are more susceptible to 
locks when the two-way side lobe level is less than 25 dB. In 
practically all cases, having a two-way side lobe level slightly 
greater than 45 dB would eliminate the possibility of a side lobe 
lock at altitudes of 4000 ft and above (for the radar parameters 
used) . 
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In conclusion: 


1) With regard to the probability of maintaining a radar 
lock, there is no significant difference between the 
modulation techniques. 

2) There is no significant difference between the modula- 
tion techniques, with respect to side-lobe power 
levels, for eliminating side-lobe locking conditions 
to an altitude of 4000 ft. 

3) Excluding the consideration of time diversity, there 
is no significant difference between the modulation 
techniques with regard to the isolation or side-lobe 
level required to prevent cross-lobe locking to an 
altitude of 4000 ft. 

4) A computer run was made with the time-correlated Monte 
Carlo parachute program using the parameters of a 

J system and the corresponding isolation determined 
o 

here (that would exclude cross-locks to an altitude 
of 4000 ft). No cross-lock conditions existed during 
the 300 parachute trajectories that were run. These 
results are shown in figures 26 and 27. 

C. Radar System Comparison 

Before deciding which modulation method to mechanize in the 
system, we first need to consider the characteristics of each 
method. This will allow us to select a method that will not be 
too complex to implement. The ICW system uses a servoed pulse 
repetition frequency (PRF) to avoid altitude holes (which would 
otherwise occur when the PRF equals the Doppler frequency) and 
maintain the proper 50% duty cycle. To maintain the 50% duty 
cycle, the slant range must be measured along the velocity beam. 
But because there is a high probability that a beam may unlock 
as a result of large attitude excursions, and an uncertainty as 
to which beam will do so, the slant range on each of the four 
velocity beams must be measured. Although this minimizes the 
need for a LARA, it increases the complexity of the four range 
channels that must be included in the basic velocity-measuring 
system to maintain the theoretical efficiency of the system. 


The reference system for the J x modulation is the APN/187 
FM/cw radar system. This particular system circumvents the need 
for frequency diversity by using time diversity, in which the 
respective velocity beams are sampled at a 24-Hz scan rate. It 
operates on a single frequency, which permits the formation of 
a multiple-beam pattern in which each beam can be sampled at the 
scan rate. The RF separation frequency from the carrier will 
also allow us to use an interleaved aperture array, in which the 
receiver and transmitter array elements are interlaced. As shown 
in Appendix D, this reduces the original 60-dB isolation between 
the receiver and transmitter antennas significantly (to between 
18 and 22 dB) . 

But because the first-order Bessel function (J x ) does not 
have the isolation of a higher-order Bessel function (J 2 or 
J 3 ) , the use of an inter leaved-antenna array becomes question- 
able. In addition, the effect of the interlaced components is 
to increase the side-lobe levels of the multiple beams about the 
center of the aperture. As the reflectivity curves reach a maxi- 
mum value about the vertical or nadir, this particular antenna 
system is more susceptible to side-lobe lockups than one that 
uses a side-by-side aperture array. 

The time-diversity system has an additional advantage: it is 

the lightest. Its disadvantages -- of using single channel 
(which creates operational as well as reliability problems) and 
its high probability of having a beam unlock, are overcome by 
using the additional RF beam to increase the probability of having 
three beams locked. Had it been only a single-channel system, 
however, the high probability of having a single beam unlock would 
cause there to be a high probability that the total system (four 
beams) would unlock . 

For application to the Mars Lander, we recommend that indi- 
vidual receiver channels be used to ensure a high probability of 
locking three beams. Due to the previously-mentioned consequences 
of having the inter leaved-antenna array operating at a single 
frequency, a time-diversity system has little or no advantage 
over a frequency-diversity system for the Mars Lander mission. 

The J q system has a weight disadvantage, because an addi- 
tional quadrature preamplifier must be provided to maintain sign 
sense; however, this disadvantage may be offset somewhat by using 
solid-state transmitters (these might be susceptible to either 
AM or FM bandwidth and transient requirements in the ICW or FM/cw 
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systems). The high-pass filter in the preamplifier, which atten- 
uate diode noise in the mixer and spectral impurities in the trans 
mitter, makes this system more susceptible to beam unlocks as the 
velocities decrease toward the zero Doppler region. However, this 
same preamplifier rolloff characteristic tends to minimize side- 
lobe and cross-lobe lockups, in that the unwanted signals from 
the side lobes occur in the proximity of the zero Doppler region, 
where the reflectivity curves peak at the nadir. 

Due to the criticality of providing altitude information from 
4500 ft to engine cutoff (10 ft), we studied methods of using the 
TDLR to provide a backup range measurement. As previously dis- 
cussed, the ICW radar system utilizes a range measurement at the 
upper altitudes to properly servo the PRF; but as the altitude 
decreases to about 500 ft, the transmitted pulse eclipses the 
received pulse and makes this method of modulation inoperative. 

To compensate for this effect, the ICW modulation is augmented by 
an FM/cw modulation whenever the altitude is below 500 ft. 


For this latter case, and for the Ji FM/cw modulated system, 
the range is generally determined by measuring the phase differ 
ence between the Bessel sideband of the transmitted and received 
signals . 


The S/N ratio for the FM/cw system is 


S 

N 


FM/cw 


P G A 2 a J 2 (M) 
w o n 

(4nR) 2 L KT B t NF cos 0 


where 

M = 2m sin ( 2:f “ 

\ m C 

and the phase or range accuracy can be expressed as 



where 

C = (989 feet/2*) (1 MHz/n f J 


(59) 


(60) 


(61) 


(62) 
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Upon substituting the S/N term into equation (61) , we obtain 


« 



// irA (LKTB~ NF cos 
|4£C\ t 

\ A A p • Go 

' 7 V w o 


0 


(63) 


where n refers to the order of the Bessel sideband. From equa- 
tion (63) we can see that, as the argument M of the Bessel 
function varies with respect to range, the one-sigma range error 
AR will also vary as a function of range unless one of two 

conditions are met: either the range measurement interval for 

R must be selected in such a way that the slope of the decreasing 

range is the inverse of the slope of the Bessel function J (M) 

n 

for the selected argument, or the S/N ratio must be very large 
(S/N » 1). In the latter case, the large S/N ratio extends 
through the receiver to the final phase measurement. 

The other well-known technique for measuring range involves 
generating a linear sawtooth waveform, of the form 

E(t) = A(t) cos t + %it 2 j (64) 


where 



Although the range is measured using a phase -measuring device, 
the frequency spectrum of the sawtooth waveform includes several 
Bessel sidebands, which constitute a rectangular frequency spec- 
trum. The number of significant sidebands f ^ ' is the ratio 

of the peak frequency deviation AF to the repetition frequency 
f ^ , and the contribution beyond this limit is insignificant. 

For the upper sidebands, the spectrum function can be expressed 
as 


F (M) 


u 


00 


E e n <« 

n=0 


(65) 
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which shows that the S/N expression for the linear sawtooth wave 
is independent of variations in the altitude. 

Thus the range may be measured either using a single Bessel 
sideband on each beam of the velocity sensor or a single-beam 
linear FM/cw modulation; however, the linear FM/ cw modulation 
technique is essentially independent of variations in the slant 
range . 

Another consideration for the Mars Lander mission is the sus- 
ceptibility of the radar system to surface effects. Since the 
exact characteristics of the local terrain are unknown, the nominal 
and "worst-case" surface characteristics must be estimated and 
incorporated into surface reflectivity curves. Three reflectivity 
curves were used in the simulation: an upper curve, which was 

determined from the specular component; a mean curve, which repre- 
sented an unsloping, diffuse surface with no surface slope; and 
a lower curve, which represented the contribution from a sloping, 
diffuse surface. The probability density function, which describes 
the percentage of specular and diffuse surface components within 
a beamwidth-illuminated area, can then be determined from a knowl- 
edge of the Martian local terrain or an assumed reference terrain. 

A brief description of the reflectivity curves is given in Ap- 
pendix E. 
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As the altitude decreases and the diagonal distance across 
the beamwidth-illuminated area becomes comparable to the surface 
decorrelation distance of a rough terrain, amplitude variations 
occur in the Doppler time video (see ref. 10). In addition, 
exceedingly large objects, whose effective radar cross-sections 
are large in comparison to other objects within the beamwidth, 
can produce large chirp signals that have different slopes than 
those from objects at the centroid of a homogeneous illuminated 
area. The above conditions tend to produce either excessive 
fading or strong signals that cause the tracker to be pulled off 
the Doppler centroid for short periods of time. Although this 
effect is common to all forms of modulation, we observed that 
not all of the reference systems had provisions to minimize this 
effect . 

Because there is a short time constant in the TDLR tracker, 
so that it can perform over rapidly-varying trajectories, the 
tracker must have a memory circuit with a suitable time constant 
to prevent the recycling that might result from these instan- 
taneous interruptions . 



VII. terminal descent system modeling 


A, General Description 

The 6D0F terminal descent computer program that generated the 
results shown in this report has been used by Martin Marietta since 
about June, 1967. It has been under continual evolution as the 
Mars Lander mission became better defined and as the mathematical 
modeling of components such as the radar and the propulsion system 
became more detailed. The latest version of this program is being 
i“ nted as . a P art of thi s contract and shall be referred to as 

unBTPAM^-'^ flOW charts > module diagrams, symbol definitions, 
and FORTRAN listings are documented in reference 1. This chapter 
summarizes the program, the mathematical models used, and the op- 
erational features. Some of this material is repeated in reference 


The functional part of the program, which simulates vehicle 
dynamics and components, is made up of building blocks called 
modules . Though a module is actually a subroutine in the general 
programing sense, in this program the term is used to designate 
these major building blocks to separate them from the executive 
subroutines and other miscellaneous subroutines. There are five 
groups of modules: 

Aerodynamics: A1 thru A5 

Control or computer: Cl thru CIO 

Dynamics: D1 thru D5 

Geodetic: G1 thru G6 

Sensors: SI thru S10 

The user may place any FORTRAN models or operations he chooses 
in any of these modules and designate the modules to be executed 
in any run by using input data cards. 

The executive structure of the program consists of the main 
program and the subroutines required to set up and terminate runs, 
read input data, write output data, set up plots, and perform 
other housekeeping functions. The user very rarely needs to get 
into the details of the executive routines. Occasionally, he may 
need to use certain of the executive subroutines to switch from 
one phase of the run to another, such as switching from the para- 

chute-descent phase to vernier engine ignition or terminal-descent 
staging . 
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Figure 38 shows a flow block diagram of the program as it is 
used during terminal descent. Starting with G1 , the modules are 
called in the order shown to compute the derivatives of the dif- 
ferential equations and the other related variables. Then the 
executive calls the integrator (FOMS/AMRK), CNTR, OUPT , etc, 
as shown. The integrator can use either a simple first-order 
integrator with derivative averaging (FOMS) or fourth-order 
Adams-Moul ton or Runge-Kutta integrators (AMRK) by specifying 
the form during his input. He may also choose between generating 
random numbers having either a uniform or Gaussian probability 
density function for use in any module (e.g., for radar noise and 
output quantization) . 

If he desires, he may include lander aerodynamics by calling 
Modules G3 , Al, and A2 into the simulation between G1 and G5 . 

The flexibility of the program is demonstrated by noting that the 
use of Modules G1 , G3 , Al, A2, G5 , Dl, and D2 and appropriate in- 
put data will yield a 6D0F entry vehicle simulation. This has 
been done several times for Mars Lander or Venus missions, with 
start-up and checkout times of less than two weeks. 

The mathematical modeling used in Modules A, G, and D is de- 
scribed in the next chapter and in reference 1. The next section 
describes the modeling for the cw linear FM ranging-type radar. 

The modeling associated with Modules C3 , C4, and C5 has been 
covered in Chapter III, with the exception of the propulsion sys- 
tem portion of C5, which is presented in Section C of this chapter. 
In general, the modeling duplicates the system block diagram shown 
in figure 4 and is based on the following assumptions: 

1) Inertial sensor dynamics and quantization effects are 
neglected . 

2) Accelerometer measurement errors and gyro drift errors 
are constants (g-sensitive effects are neglected) , 

3) The lag portion of the lead-lag network is too high 

in frequency to be simulated economically in a digital 
simulation, so only the lead portion is included with 
the valve dynamics. 

4) All filters in the system are first-order. 

5) The valve short-term transient dynamics are assumed 
to be second order, and there are no nonlinearities, 
except for position- and rate- limiting . 

6) Time lags in the engine and feed system are too short 
to be simulated. 

7) Engine misalignments and unbalanced aerodynamics are 
ignored . 


ordinate U w Radar 

nversion I Radar I ™ filters 
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Figure 38,- Flow Block Diagram of the M0D6MV Computer Program as Used During Terminal Descent 
















B . Radar Model 


The approach used in simulating the MOD-LM radar was that of 
a black box or functional simulation, as opposed to a one-to-one 
type of modeling. The following will describe the math models 
used in the simulation; as a guide to this discussion, a block 
diagram of the simulation is shown in figure 39 (the following 
text will refer to block A, B, etc, from this figure). 

Figure 39 is representative of the flow of calculations car- 
ried out on a per beam basis within 6D0F simulation program. Fig- 
ure 40 indicates the relative beam geometry. Note that the range 
beam appears in the center of the beam pattern; this is a recom- 
mended change to the system, and the LM radar has been evaluated 
with the range beam in this position. 


1. Doppler and range frequency calculation (Block A) .- The in- 
beam Doppler frequency and range frequency are calculated by: 

(1) transforming the inertial velocity of the vehicle into body- 
axis velocities; (2) transforming the body-axis velocity compo- 
nents into antenna-beam velocity components; (3) converting the 
velocity components along the beam into Doppler frequencies; (4) 
converting the inertial position of the vehicle into the position 
in the vehicle-centered coordinate system; and (5) transforming 
the position in the body-axis coordinate system into a component 
along Range Beam 4. 


a. In-beam velocities: The in-beam velocities are obtained 

by solving equations (69) thru (75). 


V B = A V x (69) 


V B1 

= 4 B1 ' 

' V B 

(70) 

V B2 

= ^B2 

' ^B 

(71) 

V B3 

= ^3 ‘ 

— > 

• V B 

(72) 

V B5 

= ^B5 ' 

t > 

(73) 
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Figure 39,- Flow Block Diagram of the Calculations for Each Beam, Modified LM Radar Simulation 











Note : 1. Beam 4 is the range beam. Con- 

figuration angles of Beams 2 and 
3 are nominally identical to those 
of Beam 1, Beam 5 only transmits. 

2. (; = 20° 22 ' 48", A = 13° 59' 22", 

0 — 24° 33 1 , \|r = 14° 53', and 

7 = 19° 45'. 


Figure 40.- Definition of Beam Pointing Angles with Respect to Vehicle 
Body Coordinates, Modified LM Radar Simulation 



where A is the inertial- to body-coordinate-system transforma- 

— > 

mation matrix; is the vector expressing the inertial velocity 


of the vehicle; 


B 


is the velocity of the vehicle expressed in 


body coordinates; \i 


Bl 3 


X B2 3 ^B3 


and are the unit vectors 

Jdj 


along Beams 1, 2, 3, and 5 in the body coordinate system; and 


V 


Bl* 


V B2’ V B3’ 


and are the components of the vehicle's 


velocity along Beams 1, 2, 3, and 5. 
b. Range: For ranges, where 


R = D/cos 6 , 
4 

(74) 

^14 = ^ ®4 ^B5 

(75) 

c°s e 4 = n I4 • H N 

(76) 

° = x . ; N 

(77) 


— 'p 

and A is the body- to inertial-coordinate-system transforma- 
tion matrix; is the beam 4- to body-coordinate- system trans- 
formation matrix; is the unit vector in body coordinates 

along the Beam 4 pointing direction; p. . is the unit vector nor- 
mal to the terrain; is the Beam 4 pointing direction in in- 

ertial coordinates; 6 is the incidence angle for Beam 4; X 

is the vector for the position of the vehicle; D is the perpen- 
dicular distance from the vehicle to the ground; R is the slant 
range of the range beam (Beam 4) . 


c. Doppler frequencies: The simulation is based on assuming 

an infinite ground plane, which can slope in a defined direction. 
Then equation (78) is used to obtain the Doppler frequencies for 
the velocity beams. 


f 


d 


2 V Bi 


(78) 
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where V is the velocity of the beam and A is the wavelength 
Bi 

of the beam. 

d. Range beam frequency: Equation (79) is used to obtain 

the range beam frequency. 


f 


R 


KR + 



(79) 


where K is a scale factor corresponding to the frequency modu- 
lation of the range transmitter. 

2. Signal characteristics (Block B) .- This description pertains 
to Block B of figure 39; the calculations are in regard to the 
signal spectrum characteristics, the bandwidth, the power level, 
the proportion of received power in the tracker step low pass fil- 
ter ( SLPF) , and the power level required by the threshold circuitry 
(as determined by the tracker output) . 

a. Incidence angle: The incidence angle of each beam is cal- 

culated in the following manner: 


li = A B. m , where j = 1, 2, 3, 4, 5 
“ “J B J 


cos 0. = Hjj • P N 


( 80 ) 

( 81 ) 


where the subscripts refer to the calculation being made for each 
of the radar antenna beams. 


b. Radar cross section: The radar cross-section as seen by 

each beam in illuminating the terrain is calculated for use 

in the radar range equation. The reflectivity model used is based 
on the Muhleman equation (see ref. 11), which is 


( n = U k cos 9 

V (sin 9 +0: cos 9 ) 3 


( 82 ) 


where a and k depend on the wavelength, tj is the ratio of 
the surface reflectivity and that of a perfectly-reflecting sphere; 
and 1 is the incidence angle of a beam. 
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Three reflectivity models from reference 11 are shown in 
figure 41. Curve 2 was used for all the M0D6MV simulation runs 
described in this report (i.e., a = 0.035, K = 0.037, and 
T} = 0.063). Radar evaluation is extremely dependent on the re- 
flectivity model that is used. This is brought out in Chapters 
IV and VI. It is very desirable to further investigate this area 

c. Radar range equation: The radar range equation that is 

used in this simulation to calculate the power received at the 

antenna P is 
K 


P = o to 

R (4jt) 2 R 2 cos 0 

where P fc is the transmitter power and G q is the two-way an- 
tenna gain. 



Angle from vertical, deg 


Figure 41.- Reflectivity Models Proposed for the 
Modified IM Radar Simulation 


(83) 
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The parameters used in the simulation representing the 
MOD-LM radar are given in table 16. These parameters were sup- 
plied by Robert Harrington of the Ryan Aeronautical Co., San Diego, 
California . 


TABLE 16.- PARAMETERS USED IN SIMULATION 
REPRESENTING THE MOD-LM RADAR 


Parameters 

Range 

Velocity 

Transmitter power 

-14.1 dBw 

-13.0 dBw 

Waveguide loss 

-1.8 

-1.8 

Pattern parameters 

-6.0 

-6.0 

Antenna gain 

27.2 

28.2 

Wavelength 

-19.8 

-20.6 


d. Signal power: We need to determine the signal power that 
is seen by the tracker so that we can compare it to the tracker 
threshold. The spectrum seen by the tracker may be modified some- 
what, due to lag of the tracker output (due to high tracking rates, 
the SLPFs may straddle the incoming Doppler spectrum as shown in 
figure 42) . 



Figure 42.- Spectrum Overlap in the SLPFs Due to 
High Tracking Rates 

When this situation is detected in the simulation, the 
received power calculated from equation (83) is modified as fol- 
lows : 


RT 


= P + 10 log 


ABW 

BW„ 


(84) 
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where ABW is the portion of the signal spectrum that appears in 
the tracker SLPF bandwidth when both bandwidths are idealized as 
shown in figure 43, and BWg is the Doppler signal spectrum 3-dB 
bandwidth. ° 


Doppler 



SLPF response 


Figure 43.- Idealized Spectrums for the SLPF Bandwidth 

e. Signal bandwidth: The Doppler signal spectrum bandwidth 

(see fig. 43) is calculated as follows: 


Af 

d A 


V B | sin 7 


(85) 


where A 7 is the average 3-dB beamwidth of the antenna, |v 

is the total velocity of the vehicle, and 7 is the angle between 
the centerline of the beam and the velocity vector, and the range 
signal spectrum bandwidth is defined by: 


2V A 7 

Af r = (SF) R A 7 tan 9 + — ^ ( 86 ) 

(The derivation is found in Appendix F.) 

The spectrum bandwidth is assumed to be proportional to 
the component of velocity perpendicular to the beam-pointing di- 
rection, as modified by the 3 -dB beamwidth of the antenna. The 
spectrum is assumed to have a Gaussian distribution. 
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f. Tracker threshold level: The simulation made use of the 

receiver sensitivity curves supplied by Robert Harrington of the 
Ryan Aeronautical Co. These curves are shown in figures 44a thru 
44d to determine the adequacy of the ratio required for tracking 
(referred to the antenna) . Figure 44a shows the sensitivity of 
the velocity receiver when the LM radar operates in the wideband 
mode. This curve can be described as having a Doppler spectrum 
centered at 5000 Hz with a bandwidth of 390 Hz and power level 
(as calculated by the radar range equation and modified by the 
amount of spectrum overlap in the SLPF) exceeding -139 dBw --and 
satisfies the track criterion. Figure 44b shows the sensitivity 
of the velocity receiver when the LM and Modified LM radars oper- 
ate in the low mode. Figure 44c shows the sensitivity of the 
range receiver. The curves shown in figures 44a thru 44c were 
originally derived to satisfy the LM radar test requirement. 

Figure 44d shows the sensitivity of the velocity receiver for 
wide-band operation when the bandwidth of the SLPFs in the Mod-LM 
is changed from 2800 to 600 Hz. The curves shown in figures 44a 
thru 44d were curve-fit with straight line segments for use in 
the digital simulation. Each of the curves is then represented 
by equations of the form: 

P TH = g ( f d t ’ f Rt) 

where P is the threshold tracking criterion referenced to the 
TH 

antenna, in dBw, and and are the tracker output fre- 

quencies in the Doppler and range channels, respectively. 

3. Mode switch (Block C) .- Mode switching is based on range and 
'velocity tracker outputs, and affects the SLPF bandwidths, tracker 
search limits, search rates, range scale factor, and receiver 
sensitivity. Table 17 describes the mode switching criteria. 

The criteria numbered 1, 2, 3, and 4 in table 17 apply for the 
LM-radar mode switching criteria established by Robert Harrington 
of Ryan Aeronautical Co. (ref. A13) , and the criteria numbered 
1-a , 2-a , 3-a, and 4-a apply for the suggested mode criteria 
change from 2500 to 1250 ft. The latter are presently being used 
in the simulation. 
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Signal level at receiving antenna, dBw 







(c) Range Channel Acquisition Sensitivity 





.000 10 000 100 000 

Frequency, Hz 

(d) Wideband Mode, Modified SLPF Bandwidth = 600 Hz 


Figure 44*- Concluded 




TABLE 17.- MODE SWITCHING CRITERIA 


Criteria 

Mode 

1. Range < the frequency equivalent 

of 2500 ft 

1-a. Range < the frequency equivalent 
of 1250 ft 

Low mode throughout the 
system 

2. Range > the frequency equivalent 

of 2500 ft 

2-a. Range > the frequency equivalent 
of 1250 ft 

High mode except the ve- 
locity in any individual 
velocity tracker between 
the equivalent frequency 
limits of ±65.5 fps will 
switch the SLPF band- 
widths to the low mode 
equivalent to that chan- 
nel 

3. If low mode has been achieved, 

then: 

Range > the frequency equivalent 
of 3100 ft 

3-a, If low mode has been achieved, 
then: 

Range > the frequency equivalent 
of 1850 ft 

High mode applies as de- 
scribed above 

4. Mode inhibit if the range beam 

or either of its compensating 
beams are in search mode 

4-a . Same as 4 . 

Status quo 


4. Noise generator (Block D) The statistical characteristics 
of the tracker output are obtained by adding the appropriate jit 
ter or noise to the Doppler frequency f^ and the range fre- 
quency f^. The ms of the frequency jitter c and the time 

constant of the noise t are functions of the Doppler spectrum 
bandwidth (ref. 12). Therefore 











It was assumed that the noise followed a Gaussian distribution. 
Each velocity channel as well as the range channel has its own 
noise generator, which produces variations on the respective sig- 
nals according to the parameters defined by equations (88) and 
(89). 

5 Bias (Block E) Although the simulation has the capability 
to add biases to the signals (see fig. 39), no runs were made 
with biases, because the outcome is predictable and is not a 
problem in any sense . 

6. Tracker (Block F) .- The closed-loop transfer function of the 
tracker is represented by a first-order lag, i.e., 

F(S) cl = 1/(1 + tS) (90) 

where 0.05 < t < 0.07. The closed-loop time constant t of the 
tracker is specified to be within these limits. The value used 
in the simulations to date is t = 0.06. 

7, Threshold and search control ( Block G) .- The threshold cir- 
cuit indicated in this block merely takes the threshold-power- 
level criterion indicated from the sensitivity curves [see eq . 

(87)] and compares it with the received power indicated in the 
tracker, [see eq . (84)]. Then, the functions described in table 
17 are carried out depending on the results of the comparison 
just described and the current state of the tracker. 

Additional logic is provided to detect the presence of cross- 
coupling between the velocity beams. This logic operates by 
checking for cross-coupled energy for the three possible cases 
for each beam. If any cross-coupled signal is above the threshold, 
the signal is processed through the tracker bandwidth, the level 
of each signal in the bandwidth is compared, and the tracker is 
biased to follow the strongest signal present. 

Cross-locking may occur for various combinations of surface 
conditions and vehicle attitudes, or may be forced to occur by 
selectively depressing the signal level in a beam for a fixed 
time and allowing the cross-coupling logic to operate. 

The 0.5-sec lockout time mentioned in Tracker State 2 of table 
17 refers to the time allowed the tracker to observe a signal 
below the threshold power level before it begins the search mode. 

The 0.1 sec mentioned in Tracker State 4 of table 18 is the time 
that the received signal must remain above the threshold (during 
the search mode) before normal tracking begins. 


TABLE 18.- FUNCTIONS TRIGGERED AS A RESULT OF THRESHOLD CRITERIA 


Tracker state 

Threshold criteria 

Function 

1 . Tracking 

p > p a 

RT TH 

Continue tracking 

2 . Tracking 

P < P 

RT TH 

Commence 0.5 sec lockout. 

Tf P RT > P TH - 3 dBw, the 

tracker is held at the last 
value it had during the lock- 
out time. If P RT < P TH - 

3 dBw, the tracker continues 
tracking during the lockout 
time # 

3 . Search 

P < P 

RT TH 

Continue search 

4. Search 

P > P 

RT TH 

Commence 0.1 sec proof 
time . 

P„ m = received 
R T 

antenna ; 

power seen by the tracker as referred to the 

P TH = receiver threshold sensitivity as referred to the an- 
tenna . 


a. Search control: The search mode occurs in the simulation 

on an individual-channel basis, according to the radar mode. The 
search limits and rates described in table 19 are for the exist- 
ing LM radar; those in table 20 are for the Modified LM radar. 

It takes 0.2 sec for all beams, regardless of radar mode, to drive 
the tracker VCO from the low limit to the high limit of the search 
range . 

b. Channel status: The status of output data is flagged in 

the following manner: if all velocity channels are tracking, a 

discrete is initiated to indicate this; if the range beam and its 
velocity-compensating beams are tracking, a discrete is also avail- 
able to indicate this. These two discretes now serve as LM radar 
data flags. A suggested change, which was incorporated in the 
simulation, is to use a discrete for each individual channel, to 
indicate tracker locks and searches. 
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TABLE 19.- LM RADAR EXISTING SEARCH RATES AND LIMITS 
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8, Data output (Block H) . - The velocity tracker outputs are con- 
verted to body-axis velocities by the following relationships: 


( 91 ) 


V 

u 


V 

v 


V 

w 


A 

4 cos A cos £ 

4 sin A (Dl ' 

\ 

4 cos A sin t 


(Di + D 3 ) 
D 2 ) 


(D 2 - d 3 ) 


(92) 

(93) 


where A and £ are the beam configuration angles defined in 

figure 40, Di > D 2 , and D 3 are the Doppler frequencies of 

Beams 1, 2, and 3, respectively, and V , V , and V are the 

u v w 

body-axis velocities described in figure 40, 


As shown in equation (94), the velocity trackers of Beams 1 
and 3 are used to correct the velocity component of the Beam 4 
range-tracker output. 


R = 




(D x + D 2 ) 




(94) 


where f R is the range tracker output, is the range scale 

factor (2.32 for high mode, 11,6 for low mode), R is the range 
along the u axis (see fig. 40), Ai is the wavelength of the 
range transmitter, and A 2 is the wavelength of the velocity 
transmitter . 


For radar-prime-mode simulations, the radar outputs are fil- 
tered by analog filters with first-order time constants. For 
radar-aided, inertial-navigator simulations, a digitized output 
is used, which has a specified count time. A quantization error 
of ±1 Doppler cycle is added by using a uniform, random-number 
generator . 
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C. Propulsion System Model 

The model of the propulsion system consists of equations that 
describe the throttle valve dynamics, the engine thrust, the spe- 
cific impulse, the flow rate, and the propellant-tank pressure 
blowdown effect. 

Figure 45 shows a block diagram of this portion of the system, 
beginning with the mixing matrix inputs, e Q , e^, and . Each 

input is filtered by a lead-lag network for inner-loop control 
compensation. It is difficult to simulate the pole frequency, P, 
because it. is about 300 rad/ sec, so in the digital model it is 
neglected. Since the mixing matrix is a linear operation, the 
zero (S+Z) is included in the linear portion of the valve model. 
The valve is essentially a second-order dynamic device, which, 
with the zero from the lead-lag network, is described by 



Furthermore, when X is on either the upper or lower limit and 

•L • 

is trying to drive harder into the limit, is set to zero. 

Valves 2 and 3 are described by the same equations, with the ob- 
vious change of subscript. 
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Valve Engine 



















The propellant tank pressure blowdown is assumed to occur 
adiabat ically , i.e., without the addition of heat, which is a 
conservative assumption. The initial pressurant gas volume is 
calculated from equation (100) to give the right blowdown pres- 
sure ratio when the input fuel is 100% used. 

P I- p o ( 1 - 0 + fSti; - c m )' 5 < 10 °> 

where 

C pT1 = (BLD) 1/& - 1 (101) 

and BLD is the blowdown ratio, Fug1 t is the initial fuel load 

in lb, 6 is the ratio of the specific heats, and P is the 
initial tank pressure in psia. 

The feed system is represented by steady-state pressure drop 
equations that seek to linearize the actual flowrate as closely 
as possible to the actual position of the valve. This is done by 
calculating a valve flow resistance which, with the other pressure 
drops in the feed system, gives the desired flow for the average 
pressure in the tank. The actual flowrate, however, is a linear 
function of valve postion only at this pressure and deviates some- 
what at either higher or lower tank pressures. 

In the simulation, we only approximated the steady- state , 
fluid-flow performance of the feed system because we did not 
account for the small changes in the thrust coefficient that 
occur at low flowrates. The equations used in the simulation 
are : 

a. Desired flowrate: 

= W B * X l (102) 

b. Flow resistance: 

'2 

(103) 
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where : 
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(108) 


where W is the design flowrate in lb/sec, is the blow- 

JdD 

down tank pressure in psia about which to linearize the flow 


resistance, D 


PVM 


is the valve pressure drop in psi that occurs 


at the maximum flowrate, W is the maximum flowrate in lb/sec, 

max 5 

? CM s max i mum chamber pressure in psia, CD is the injector 

discharge coefficient, and Cp is the fuel density in slug/ft 3 . 

c, Actual flowrate: 


w i = 


W C1 + 


n — 

,W C1 +2 

w„ l 
L C2 

(»c 3 ) 

1 + «C3 1 

(k 2 )1 


(PT) 


(109) 
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Percent flowrate: 
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( 110 ) 


max 
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Specific impulse : 


Igp = Table lookup for I^p as a function of 

(This table lookup is based on test data that cover the 
range of flowrates at which the engine can operate.) 


f. Engine thrust: 


= w i M 


If desired, a constant thrust or regulated propulsion sys- 
tem can be used in the simulation by setting the blowdown ratio, 
BLD, less than one. 
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VIII. PARACHUTE DYNAMIC S 


A. General Discussion 

The 6D0F equations of motion for a parachute -Lander combina- 
tion in the presence of winds and gusts are described below. In 
this chapter, the results obtained with this simulation are com- 
pared with those from an earlier two-dimensional program (ref. 13). 
The simulation of the parachute descent begins with the parachute 
open and descending under trim conditions. When the Lander reaches 
vernier ignition altitude, the vernier engines are allowed to warm 
up for a few sec before the parachute is released. When the para- 
chute is released, the parachute equations cease being computed, 
and the simulation focuses on the performance of the Lander cap- 
sule and its radar system. 

The six basic assumptions inherent in this simulation are 
listed below: 

1) The canopy and shroud lines act as a rigid body; 

2) The canopy is symmetrical and hemispherical; 

3) The elongation of the riser is proportional to the 
load; 

4) The riser is massless and transmits tension and tor- 
sion only; 

5) The capsule attachment harness acts like a rigid 
truss as long as the riser pull angle, -X , is less 

than half the apex angle of the truss (see fig. 46); 

6) The body-axis system corresponds to the principal 
axes . 



Figure 46.- Configuration of the Capsule 
Attachment Harness 
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Some of the main features of the simulation are I 

1) Both the parachute and the capsule have full 6D0F; 

2) The aerodynamic forces and moments on the parachute 
and capsule can be expressed as functions of the 
angle -of -attack and/or the Mach number; 

3) Horizontal and vertical wind profiles can be expressed 
as functions of the altitude and/or random wind gusts; 

4) The apparent and enclosed masses of the parachute are 
included; 

5) The simulation can be run with or without a parachute; 

6) Coriolis terms are available; 

7) The size of the parachute can be expressed as a func- 
tion of time in order to simulate the reefing of the 
canopy; 

8) The gravity is expressed as a function of the alti- 
tude ; 

9) The atmospheric density and the speed of sound are 
functions of the altitude; 

10) Spherical- or flat-planet options are incorporated. 


B. Development of Model 

Coordi nate system .- A plane t -fixed , tangent-plane coordinate 
system ( x Tp » Y pp , Z xp) initia ^Y oriented with respect to 

the rotating planet, as shown in figure 47. Note that the longi- 
tude of the planet is not defined. 



Figure 47.- Orientation of the Tangent-Plane 
Coordinate System 



Figure 48 shows that an Euler angle transformation is required 
to express the tangent-plane coordinates of each body as body-axis 


coordinates 


(V 


V 


| about the axisj , 


, where is the first rotation 
is the second rotation (about the Y 1 


axis), and cp is the third rotation (about the X 11 axis) 



Figure 48.- Euler Angle Transformations Needed to Express Tangent-Plane 
Coordinates as Body-Axis Coordinates 


2. Equations of motion .- The forces on each body are computed 
in the body-axis coordinate system of each body. A set of body- 
axis acceleration equations [eq. (112) thru (114)] is then re- 
solved into a set of planet -f ixed , tangent -plane equations, which 
can then be integrated to obtain the velocities and positions for 
each body. 
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( 112 ) 

(113) 

(114) 


(115) 

(116) 

(117) 

where: F , F , and F are the aerodynamic forces; 

*ba y ba z ba 

T , T , and T are the thrust forces; F , F^, and 

^BA Y BA Z BA 

are the coupling forces transmitted to the capsule or chute 

by the riser (expressed in the body-axis system); L , M BA’ 

and N- 3 - are the aerodynamic moments; L , M , and N are 
BA i 1 X 

the thrust moments; L , M^ , and N z are the coupling moments 

produced by the riser force about the roll, pitch, and yaw axes; 
and A, B, and C are the principal moments of inertia about 
the roll, pitch, and yaw axes. 

3. Determination of coupling force .- After determining the center- 
of -gravity positions and the Euler angles for the capsule and the 
parachute, the body -axis locations (denoted by the subscript B) 
of the riser attach point for the parachute (P) and the riser 
attach point for the capsule (C) are computed. These locations 
are then transposed to the tangent -plane coordinate system (see 
fig. 49) in order to determine the stretched riser length, SXZI . 
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4. Forces and moments transmitted to parachute .- Figure 50 shows 
the forces and moments that are transmitted to the parachute. 
Before attempting to solve equations (123) thru (132) , point 
v Y Z must be expressed in the chute-body coordinate 

CP’ CP’ CP 

system. K , the torsional spring constant of the riser, [eq. 
(129)] is equal to zero if a swivel is used. 
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cos A = I X^,„ - X I 
CP PBA 


sin B = |Y cp - Y pBA |/l 


sin c = ( z cp Z PBA 


SXZI 

f 

SXZI 

r 

SXZI 


a. Force components 


F = T 
XP R 


( x cp ' x pba)/ 


f yp t r y cp ' y pba 


f zp t r( z cp “ z pba)/ 


SXZI 

I 

SXZI 

r 

SXZI 


b. Moments in body axis: 


L XP = K s 


“yP 


N = F 
ZP YP 


(%) 

p ( lp ) 


%-/(*- p d) 


dt 


5. Forces and moments transmitted to capsule .- Figure 51 shows 
the forces and moments that are transmitted to the capsule. Here 


Points |X, 


and (X cp , Y cp , Z cp ) must be 


. "CBA’ ^CBA ’ °CBa) 
expressed in the capsule-body coordinate system before attempting 
to solve equations (133) thru (142). 


(123) 

(124) 

(125) 

(126) 

(127) 

(128) 

(129) 

(130) 

(131) 

(132) 
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( 134 ) 


sin C 


Z - Z 
CBA PC 

SXZI 


( 135 ) 
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6. Wi nd and gusts .- Wind may be input into the problem either 
with a steady wind profile that has horizontal and vertical 

(W v ) components, or with random gusts, using the probabilistic 

capability described in reference 2. 



( 143 ) 

( 144 ) 

( 145 ) 

( 146 ) 

( 147 ) 
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Figure 52.- Orientation of the Steady-Wind -Velocity Components 
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a. Total wind: Random wind gust components are added to 

^sy’ an< ^ W SZ ^ called for. 


w = w + w 

TX W SX GX 


w = w + w 

TY SY GY 


(150) 

(151) 



(152) 


Figure 53.- Orientation of Wind-Velocity Components with Respect 
to the Air Mass 
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7. Aerodynamic forces and moments on parachute . - The aerodynamic 
forces and moments on the parachute are shown in figure 54. In 
this figure, C is the axial force coefficient, C is the 

normal force coefficient, is the pitch or yaw moment coef- 
ficient, C L is the roll moment coefficient, is the pitch 

or yaw damping coefficient, and C^p is the roll damping coef- 
ficient . 



Figure 54.- Diagram Showing the Aerodynamic Forces and Moments 
on the Parachute 
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Note that if 

SQNV = 0, then 



sin i] = cos rj = 0 

(161) 

but if SQNV + 0; 

then 
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c. Typical Module Calling Sequence 


The modules are typically called in the order shown in figure 
55. The symbols shown in parentheses are the parameters which 
are calculated in the modules. 



Figure 55.- Flow Block Diagram of the Order of Computations 
in a Typical Module Calling Sequence 
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D. Simulation Results 


The equations of motion and their programing were verified by 
comparing the results of this simulation with the results from a 
two-dimensional parachute simulation in which a high level of 
confidence exists (ref. 13). The available integration routines 
were exercised to determine whether they could be used and to 
evaluate their effectiveness. Various computing intervals and 
riser spring constants were examined to demonstrate that realistic 
values of these parameters were executable. Each of the major 
program options was checked for proper operation, and before 
terminating the run, the parachute-release logic was verified. 

The example used in the simulation was that of a 48-ft-diam- 
eter parachute decelerating a QOO-lb capsule in the minimum 
H p,S Mars atmosphere. The program began at an altitude of 16 000 

ft, shortly after the aeroshell was jettisoned. The wind profile 
used in this simulation is shown in figure 56. 


24 x 10 3 



Figure 56.- Mars Minimum H Atmosphere Wind 

P 7 S 

Profile Used in the Parachute- 
Lander 6D0F Simulations 
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Although the wind profile shown in figure 56 does not agree 
exactly with the latest wind specification given in reference 8, 
it was used in order to make the simulation compatible with the 
single-plane program, which had been executed before the latest 
specification was issued. However, before the wind profile was 
used in later combined parachute -capsule terminal descent simula- 
tions, it was updated to reflect the minimum g atmosphere 

defined in reference 8. 

1. Comparison with two-dimensional program .- The results from 
the MOD 6MV program are compared to those from the single plane 
program of reference 13 in table 21. 


TABLE 21.- COMPARISON OF THE RESULTS OBTAINED FROM THE 
M0D6MV AND SINGLE -PLANE COMPUTER PROGRAMS 3 


Parameter 

MOD 6M V program 

Single-plane 

program 

Relative velocity, fps 

165.2 

164.17 

Inertial velocity, V^, fps 

-285.4 

-285.8 

Inertial velocity, , fps 

164.6 

163.75 

Altitude, ft 

9000 

9021.6 

Flight path angle, deg 

-29.9 

-29.8 

Chute riser force, lb 

351.7 

349.5 


a The parachute check runs were terminated arbitrarily at 9000 ft, 
after passing through two wind-sheer spikes. The data shown in 
tables 21 thru 23 simply summarize the conditions at the end of 
each run. 

The response of the system is shown in figure 57a. Note that 
there is good agreement with the data obtained from the single- 
plane program. The curves in figure 57b show that the capsule 
attitude rates are very nearly equal in magnitude and have similar 
frequency patterns . 
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(a) Altitude vs capsule attitude 


Figure 57.- Parachute-Lander 6D0F Simulation Results: Wind 

Gust Response Comparison of the Wind Gust Res- 
ponse Obtained from the Single-Plane and the 
M0D6MV Parachute-Lander Simulations 
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Capsule attitude rate, deg/sec Capsule attitude rate, deg/sec 







2. Integration routines .- Two integration routines were evaluated 
in terms of their accuracy and the computing time required. Table 
22 compares the results obtained using simple, step-wise integra- 
tion (FOMS) and fourth-order, Runge-Kutta integration for 43.6 
sec of real time. 


TABLE 22.- COMPARISON OF THE RESULTS OBTAINED USING 
DIFFERENT INTEGRATION TECHNIQUES 


Parameter 

Integration 

technique 

Runge-Kutta 

FOMS 

Fixed computing interval, sec 

0.005 

0.005 

Computing time, sec 

569 

179 

Relative velocity, fps 

165.11 

165.13 

V x , fps 

-285.3 

-285.38 

V fps 

164.51 

164.54 

Altitude, ft 

8999 

9001 

Flight path angle, deg 

-29.9 

-29.9 

Riser load, lb 1 

351.6 

351.3 


This table shows that the use of the more-sophisticated Runge- 
Kutta integration routine is not justified, in view of the much 
longer computing time that it required. 

3 . Comput ing interval . - The parachute simulation described in 
reference 13 used variable computing intervals, with a lower 
limit of 0.00312 sec. To achieve good agreement with the re- 
sults of that simulation, a fixed computing interval of 0.005 
sec was used in the parachute-lander simulation. Although a 
0.005-sec computing interval may be required for other phases of 
the terminal descent and landing system, table 23 shows that a 
computing interval of 0.01 sec is adequate for the parachute 
phase. 
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TABLE 23.- COMPARISON OF THE RESULTS OBTAINED USING 
DIFFERENT COMPUTING INTERVALS 


Parameter 

Comput ing 

interval 

0.01 sec 

0.005 sec 

Computing time, sec 

119 

179 

Relative velocity, fps . 

165.217 

165.13 

v x> f P s 

-285.48 

-285.38 

v f p s 

164.61 

164.54 

Altitude, ft 

9001 

9001 

Flight path angle, deg 

-29.9 

-29.9 

Riser load, lb 

351.55 

351.3 


In general, the computing interval should be small enough to 
follow the frequency response of the system. For this particular 
problem, the period of the longitudinal oscillations in 0.33 sec 
(based on a spring constant of 1000 lb/ ft) and the period of the 
angular capsule oscillations is 2.05 sec.. If 20 computations are 
to be performed per cycle, the maximum computing interval will be 
0. 0165 sec. 

4. Spring constant .- The stiffness of the equivalent spring con- 
stant of the riser has a bearing on the computing interval re- 
quired. The stiffer the spring, the higher the longitudinal fre- 
quency, and therefore, the smaller the computing interval required. 
The checkout runs included herein were made using a spring con- 
stant of 200 lb/ ft, since little difficulty is experienced with a 
soft spring. Subsequent runs used a more realistic value of 1000 
lb/ ft to verify that the program can handle the higher frequency 
involved. 

5. Pitch and yaw case. - The motion of the system, which involves 
all the equations of motion, was evaluated by moving the wind di- 
rection 45 deg out of the pitch plane. (All the previous compari- 
son runs, however, were run with the wind in the pitch plane.) 

The wind profile and altitude time-histories in figure 58 were used 
to correlate the gust response with altitude and time. 

The capsule pitch- and yaw-rate responses to wind gusts are 
shown in figures 59a and 59b. The responses begin at four sec and 
again at 36 sec. The magnitude and frequency of these curves are 
nearly equal -- as they should be. The pitching-over of the tra- 
jectory appears as a negative bias in figure 59a. 



Wind velocity, fps 
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The pitch Euler angle time-history is shown in figure 60a. 

Except for wind spikes (which produce bumps on the curve beginning 
at four and 36 sec) the curve reflects the trajectory pitchover 
from -0.44 rad (-26.5 deg) to -1 .57 , radians (-90 deg). The yaw 
Euler angle plot (fig. 60b) has large swing after 36 sec. This 
results because the yaw Euler angle switches quardants after having 
pitched over 90 deg. To fully understand this result it is neces- 
sary to use figure 60c to trace out the ordered Euler angle rota- 
tions . 

6. Roll response .- To check out the roll coupling between the 
parachute and the capsule, we gave the capsule an initial roll 
rate of 1 rad/sec and used torsional spring constant of 10 ft-lb/ 
rad and a chute roll-damping coefficient of -0.1 rad/sec. The 
resulting capsule and parachute roll response shown in figure 61 
indicates that the chute will ultimately reduce the roll rate of 
the capsule without picking up much roll of its own. Although 
the response shown here may not be typical of parachute behavior, 
the system is responding correctly, in accordance with the equa- 
tions of motion. 

7. Combined parachute- terminal descent simulation results .- As a 
final checkout of the parachute program, a second run was made. 

This run began with the capsule at an altitude of 16 000 ft. 

After it descended through three wind gusts, the vernier engines 
were started and the parachute was staged. Then the capsule com- 
pleted the terminal descent. 

Some of the plots obtained from this run are shown in figures 
62 thru 67. The radar was initialized in the unlocked condition. 

It locked and unlocked several times due to wind gusts during the 
parachute descent, but locked and remained locked during the 
terminal descent. This can be seen in the plot of U-axis ve- 
locity shown in figure 63. The effect of the 25-deg surface slope 
can be seen in the altimeter data of figure 62. Figure 64 shows 
how the pitch rate and pitch attitude varied with time. Notice 
that the maximum rate response to the wind gusts was 22 deg/sec. 

The thrust history for Engine 1 is shown in figure 65. During 
the parachute-descent phase, this engine if off; but at T = 50 
sec, the engines are started, warmed up for two sec with the 
control loops open, the control loops are closed, and then the 
parachute is released. The thrust spikes that occur during the 
terminal descent phase result from attitude maneuvering and 
radar noise. 

The altitude velocity curve, figure 66, shows the effect of 
the wind and gusts that are plotted in figure 67. During this 
run, the lander aerodynamics were not omitted after the parachute 
was released. Thus, after the engine started the velocity de- 
creased significantly, due to the additional drag force. This ef- 
fect is normal, if we assume that the engine thrust does not change 
the lander aerodynamic coefficients. 

Martin Marietta Corporation 

Denver, Colorado, August 25, 1969. 



Capsule pitch rate, rad/sec 
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(a) Pitch rate response 


Figure 59.- Parachute -Lander 6D0F Simulation Results: Capsule 

Pitch- and Yaw -Rate Responses to Wind Gusts 
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Capsule pitch Euler angle, rad 
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(a) Pitch-Euler-Angle Response 


Figure 60.- Parachute-Lander 6D0F Simulation Results: 
Capsule Pitch- and Yaw-Euler-Angle 
Responses to Wind Gusts 
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Capsule roll Euler angle, rad 



(c) Roll -Euler -Angle Response 
Figure 60. -Concluded 
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Figure 62.- Second Parachute-Lander Simulation Results: 
Altitude vs Time 
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Figure 64.- Second Parachute-Lander 6DOF Simulation Results: 
Pitch Rate and Pitch Attitude vs Time 
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Figure 65.- Second Parachute -Lander Simulation Results: 
FI Engine Thrust vs Time 
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Figure 66.- Second Parachute -Lander 6D0F Simulation Results: 

Descent Velocity vs Altitude 
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APPENDIX A 


RADAR-AIDED INERTIAL NAVIGATION 
USING A KALMAN FILTER* 


*The work presented in this appendix is based on notes 
written by Robert N. Ingoldby, the principal investigator for 
this portion of the study. 



APPENDIX A 


A Kalman filter is one method of formulating a statistical 
estimation process to extract information about the states of a 
dynamic system from measurements of signals related to the states. 
When applied to linear systems in which Gaussian white noise 
corrupts both the system and the measurements, the theory is 
well developed and has been applied, in concept at least, to 
many problems involving guidance, control, and information proces- 
sing. In these cases, the filter gives a minimum-variance esti- 
mate . 


If the system is nonlinear, as is the usual case, then the 
linear theory can still be applied if perturbation techniques 
are used to obtain linear approximations of the nonlinearities 
in the system. 

If we make these assumptions, the dynamic equation for state 
propagation is 

Vi ‘ W( x k) + \ 

and that for the measurements is 

Z k+1 = \+l X k + ^k+l 


where 


X. is a vector of N states at time t = t, , , 

k k k+l,k 

is the state transition matrix that propagates the states at 


t, to states at 
k 


J k+1 


k+1’ 

is a vector of M 


is the dynamic modeling noise, 


measurements that is related to the 


states by the matrix 


and 


! k+l 


is the measurement noise. 


The Kalman filter makes use of the measurements and the known 
dynamic equations to generate an estimated state, X. The 
equations for the filter are represented in discrete form as: 

\+l = ®k+l,k \ + hc+l [ Z k+1 " \+l \+l,k \] 


where 



k+1 


= $ P 

k+1 ,k k 


\+l,k + Q k+1 


(Al) 

(A2) 


(A3) 

(A4) 

(A5) 
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r.* 


k+1 P k+1 ‘ \+l \+l P k+1 


(A6) 


Let us examine these equations. In equation (A3) we see that 
the new state, , is generated from the old state, , using 

the dynamic model, $^+1 k y an< * a term t ^ at includes the new 

^k+l k is w ^ at t ^ le fi^- ter 

thinks the new measurement will be, the bracketed term in equation 
(A3) is an error term, usually called the measurement residual. 

This error is multiplied by the optimum gain matrix to 

determine the correction term that must be applied to the esti- 
mated state. 


measurement, Because 


Equations (A4) thru (A6) define how the gain matrix is com- 
puted . This matrix is based on the statistics of the measure- 
ment noise, the dynamic modeling error, and the initial esti- 
mates of the state. is the covariance matrix of measure- 

ment error, 

\ ■ c ° v [\ \] 

is the covariance matrix of the dynamic modeling error 

\ ■ c ° v b <] 


(A7 ) 


(AS) 


and P, is the covariance matrix of the errors in the state 
k 

es timate 

(A9) 

Matrices and and the initial conditions of P^ can 

not be defined without a knowledge of the particular system in- 
volved . 


P. = Cov 
k 


(\ - \) (\ - \) T ' 


For the radar-aided inertial navigation problem, the states 
of interest are: the three body axis velocity components, the 

altitude, and the direction cosines relating the gravity vector 
to the body-axis coordinates. These states are described by 
the navigator equations: 
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u = a x - q . w + r • v + g . A l3 ( A10 ) 

V = A y + p . w . r . u + g . a 23 ( A11 ) 

w = A z - p . v + g . u + g • A 33 ( A12 ) 

• 

H = ” a 13 • u - A 23 • V - A 33 . W (AI3) 

a 13 = a 23 * r - A 33 . q ( A1 4 ) 

A 23 = A 33 • P - A 13 • r (A15) 

A 33 = A 13 ' q - A 23 • P (AI6) 


where u, v, and w are the body-axis velocity components 
along the X fe , Y fe , and Z axes, respectively, A , A , and 

A z are the measured accelerations in these directions, p, q, 

and r are the attitude rates about these axes, g is the 
gravitational constant assumed for the planet, H is the alti- 
tude of the vehicle with respect to the planet, and Aj 3 , A 23 , 

A 33 are direction cosines. 

If equations (A10) thru (A16) are mechanized in the flight 
computer, they can replace the state propagation matirx, 1>, 
in calculations made to update the states with time. 

If the Kalman filter is to be mechanized to estimate only the 
four .states, (u, v, w, and H) , then the measurement matrix, 

H, is the identity matrix, because the radar measurements are 
U R’ V R’ and W R from the velocity radar and from the 

altimeter. Thus, equation (A3) can be computed from the navigator 
outputs and the measurements if we are given the gain matrix, K 

For conventional Doppler inertial systems, the elements of 
K are selected as constants, usually in diagonal form. Alter- 
natively, the Kalman filter can be used to compute an optimal 
time-varying gain matrix; in this computation, the assumptions 
of linearity and Gaussian statistics are required. 
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If in equations (A10) thru (A16) only u, v, w, and H 
are to be estimated, then all the other parameters (p, q, r, 
etc) can be considered as time-varying coefficients in equations 
(A10) thru (A13) . Then the equations are linear and $ can be 
computed using standard techniques. In this analysis, we as- 
sumed that would be computed by the flight computer using 

rectangular integration, which yields 


1 /dt 

r k 

- q k 

0 

- r k 

1 /dt 

P k 

0 

q k 

" P k 

1 /dt 

0 

-Al3 k 

- A H3 k 

-A 3 3 k 

1 /dt 


(A17) 


The covariance matrix of measurement errors, R, is nondiagonal 

because u D , v_. , and w are computed from the same radar beam 
R R K 

velocities . R has the form 


C C C 0 


uu 

uv 

uw 


uv 

C 

vv 

C 

vw 

0 

*uw 

C 

vw 

C 

ww 

0 

0 

0 

0 

c 


(A18) 


R is generally time -varying , since the radar errors are a func 
tion of velocity and altitude. The numerical values used in 
equation (A18) can be obtained from error analysis of the radar 


The covariance matrix of dynamic modeling errors, Q, con- 
tains all the errors that could significantly affect those states 
which have not been modeled in propagating the state estimates 
from one point in time to the next. These errors include: 


1) Accelerometer errors; 

2) Gyro errors; 

3) Errors in the unestimated states of A 13 , A 23 , < 

A33 » 

4) Errors due to the variation of gravity; and 

5) Errors due to the effect of surface slope on alti 
tude (referenced to the landing site). 
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Over the short trajectory times for which the filter will 
be operating, the accelerometer, gyro, and gravity errors have 
negligible effect. However, the effect of errors in the un- 
modeled states A 13 , A 23j and A 33 can be significant. In 
equations (A10) thru (A16), we see that if we perturbate the 
direction cosines as 


a 13 = A 13t + AA 13 (A19) 

A 23 = A 23 t + AA 23 (A20) 

A 33 = A 33j- AA 33 (A21) 

then after rectangular integration of the equations, error terms 
exist as follows: 


Au = g • Mi 3 • dt (A22) 

Av = g • M23 • dt (A23) 

Aw = g • M 33 • dt (A24) 

AH = -(AA 13 • u + AA 23 • v + AA 33 • w) dt (A25) 

then Q becomes 

Q = Cov Y Y T (A26) 

T 

Y = [Au, Av, Aw, AH] (A27) 


To account for the surface slope effect on the altitude esti- 
mate, a term representing its mean square error can be added to 
the (4, 4) element of Q. 


Two special modules were included in the M0D6MV computer pro- 
gram to evaluate the operation of the Kalman filter. These 
modules solved equations (A10) thru (A16) and the Kalman filter 
equations during the parachute- and terminal -descent phases of 
the Mars Lander mission. The following discussion summarizes 
the results of the computer runs. 
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Q was computed using u = 300 f ps , v = 20 fps, w = 20 fps, 
3a attitude errors (10 deg) in pitch and yaw, a nominal pitch 
angle of -60 deg, and a nominal yaw angle of 0 deg. Then R 
was computed using la errors of 1.57. along each Modified-LM 
radar beam and a la altimeter error of 32 ft. The numerical 
result was 

12.2 23.0 16.4 0 

_ _ 23.0 174.0 62.0 0 

16.4 62.0 89.0 0 

000 1024.0 

If a radar beam unlocks, the filter should be set to ignore 
the data associated with that beam. In the simulation, this was 
accomplished by increasing the appropriate elements in R until 
the radar beam relocked. This caused the corresponding elements 
in K to go to very small values, thus effectively ignoring the 
affected measurement. 

Initially, P was loaded with large values. A la velocity 
error of 300 fps and a la altitude error of 150 ft were used. 
However, in all runs, the first valid radar measurement reduced 
the values in P close to those in R. 

Five computer runs were made to evaluate the operation of 
the Kalman filter. During the parachute phase, the filter was 
merely updating the state estimates every 200 msec; but during 
the terminal descent, the estimated states (still being up- 
dated) were used for vehicle guidance in the radar-aided inertial 
navigation mode. 


Run B265 

Q and R were as defined for the simulation and were con- 
stant matrices through the entire run. As the vehicle descended 
and the velocity decreased, the filter began ignoring the radar 
data. This is because after some time the terms in P and K 
become small since the filter thought it knew what the trajectory 
was; since R was computed at high values of velocity, it gave 
an unnecessar ily -conservative estimate of the measurement errors 
at low velocities. Thus, the filter ignored the data after it had 
reduced the values in P below those in R. 
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Run G439 

In this run, appropriate terms in R were made proportional 
to velocity squared and altitude squared (this assumption is 
consistent with a percentage error specification for the radars) 
and Q was held constant. This filter became unstable because 
when the side velocities went to zero, the filter thought it was 
getting perfect data. 


'Run B631 

In this run, Q was increased by a factor of 20. R was 
computed as in the previous run, except that no term was allowed 
to decrease below 0.5. In this run, the filter relied upon the 
radar data to an excessive degree, which caused the vehicle steer- 
ing to be very sensitive to radar noise. 


Run B521 

This run is the same as Run B631, except that Q was only 
four times as large as the Q of Run B265. This run produced 
fairly good results. Although there was some lateral velocity 
(5 fps) at landing, this could probably have been eliminated 
with further juggling of the computer gains. 

In conclusion, it appears that a Kalman filter can be used 
for radar -aided inertial navigation, even though its benefits 
are of questionable value when compared to the simpler naviga- 
tion approach described in Section III-A. 

The following specific conclusions can be made: 

1) The filter will not work satisfactorily unless R 

is computed as a function of the vehicle velocities. 

2) Q must be chosen carefully and must be verified 
by simulation. 

3) When filter outputs are used for closed-loop steer- 
ing during terminal descent, they can give rise to 
stability problems that result because of difficulty 
in predicting the nature of K, which is actually 

a gain incorporated into the control loops. 

4) The results obtained here show that the simpler 
navigation approach is difficult to improve on, 
particularly, if the initial radar data is used 
directly to initialize the navigator outputs. 
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DESCRIPTION OF THE TIME-CORRELATED 
MONTE CARLO RADAR LOCK PROGRAM 



APPENDIX B 


This Monte Carlo program evaluates the lock status of the ra- 
dar beam on a time-correlated basis. The program is able to simu- 
late various parachute trajectories or terminal-descent trajec- 
tories (500 flights provides a reasonable computer cost) to study 
the effect of atmospheric densities, wind velocities, wind gusts, 
and terrain slope. In this manner, the program provides a measure 
(in a probabilistic sense) of success in meeting guidance and con- 
trol requirements. The program will provide the following data: 

1) The probability of having at least three radar beams 
locked T sec before the capsule descends to engine- 
ignition altitude (4000 ft), or engine-cutoff altitude 
(10 ft). 

2) The probability of having three radar beams locked for 
at least T sec before the capsule descends to engine- 
ignition altitude or engine-cutoff altitude. 

3) The probability of having beam cross-lock T sec before 
the capsule descends to engine-ignition altitude or 
engine-cutoff altitude. 

4) The probability of having beam cross-lock for T sec 
before the capsule descends to engine-ignition altitude 
or engine-cutoff altitude. 

For the parachute-descent time-correlated Monte Carlo program 
tables have been compiled that list altitude, velocity, and flight 
path angle as functions of time for three specific atmospheric den- 
sities. At the start of each trajectory or case, the atmospheric 
density will be sampled using a probability histogram as defined in 
reference 8. An interpolation subroutine is used to determine the 
proper values of velocity, altitude, and flight path angle. In a 
similar manner, wind velocity, gusts, and gust altitude are deter- 
mined with the aid of histograms defined in reference 8. Then, 
using the tables, a trajectory is generated at specified time in- 
tervals and the beam velocities, incidence angles, etc, are com- 
puted (see Appendix C) . The radar lock status will be’computed, but 
this time the radar search function is included. As time is incre- 
mented down a trajectory, the time during which each beam is locked 
or cross-locked is also computed. As this process is repeated, say 
500 times, an overall time-of-lock probability histogram can be 
constructed, which includes the random nature of the trajectory as 
well as time-dependent effects such as gust response and radar sweep 
time. 
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The geometry computations at each time point are handled in 
the same manner as those in the static Monte Carlo program de- 
scribed in Appendix C. 

Statistics are tracked for 3-, 4-, and 5-beam radars using 
modified LM radar parameters* A flow diagram is shown in fig- 
ure Bl. 
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Figure Bl.- Flow Block Diagram of the Parachute Descent 
Time-Correlated Monte Carlo Radar Lock Com- 
puter Program 
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Figure Bl.- Continued 
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Figure Bl.- Concluded 
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DESCRIPTION OF THE STATIC 
MONTE CARLO RADAR LOCK PROGRAM 
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This Monte Carlo program works on the premise that the results 
of a large number of specific cases, each case being randomly com- 
puted, will yield the probability that a certain outcome will oc- 
cur. For instance, the result may be that 47 out of 1000 times a 
radar beam will be unlocked at altitude of 4000 ft. Figure Cl 
shows a flow chart of the program. In general, the technique used 
is to obtain random samples for surface slope, wind velocity, etc, 
and then, for that specific set of random numbers, go through the 
geometry to compute beam incidence angles and the velocity along 
each beam. Next, the radar threshold (required received power re- 
ferred to the antenna) is computed for each beam. For the LM ra- 
dar, this threshold is a function of the velocity of the beam or 
the doppler frequency in the beam. Then, using the incidence 
angle, the altitude, and a reflectivity curve (in this case, the 
one shown in fig. 23), the radar range equation is used to compute 
the received power in each beam. 

If this power is above the threshold the beam is locked; if 
the power is below the threshold, the beam is unlocked. We as- 
sume that the condition has existed long enough for the radar 
proof time to be satisfied. The equations for these calculations 
are the same as those given in reference 4 for the LM radar. 

Furthermore, for each beam, a side-lobe power level is com- 
puted by subtracting 28.9 dB from the main-lobe power. Then 
these side lobe power levels are checked against the thresholds 
of the other appropriate beams to determine whether a side-lobe 
lock potentiality exists. For instance, assume that the thres- 
hold in Beam 2 is -150 dB and that the received power in Beam 1 
is -100 dB. Then, the side-lobe power from the Beam 1 direction 
that is seen by Beam 2 would be -128.9 dB, and a potential cross- 
lock situation would exist. Further assume that the main-lobe 
power in Beam 2 is -160 dB, so that an unlock in Beam 2 is de- 
clared. The Monte Carlo program would then increment a counter 
to record that Beam 2 was unlocked at the same time that a po- 
tential cross-lock condition existed. Taking many cases, this 
would produce the joint probability function P(c,u) discussed 
in reference 14. Counters are included to record all combina- 
tions of beam unlock and cross-lock, including cross-lock with 
the transmitted fourth beam of the LM radar. 

Notice that the existence of a potential cross-lock condition 
does not mean that the radar will definitely cross-lock. The 
probability of cross-lock depends on the time required for the 
tracker to move to the side-lobe frequency, the proof time, and 
whether the geometrical conditions creating the side lobe-power 
remain there that long. This probability can be assessed by the 
time-correlated Monte Carlo program. 
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Figure Cl.- Flow Block Diagram of the Parachute-Descent Static Monte 
Carlo Radar Lock Program 
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A. Monte Carlo Program Geometry 

Figure C2 shows the geometry used in the program to calculate 
beam incidence angles and velocities. First assume the following 
coordinate system on the surface: 


Vehicle 



Surface 


Figure C2 . 

- Geometry Used in the Parachute- 
Descent Static Monte Carlo Radar- 
Lock Program to Calculate Beam 
Incidence Angles and Velocities 


The coordinates 

of the unit vector normal to the slope in 


the 1,2,3 system 

are : 



YSLOP(l) = COS (SLOPE) 

(Cl) 


YSLOP (2) =0. 

(C2) 


YSL0P(3) = -SIN (SLOPE) 

(C3) 
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The vehicle trajectory may be going up or down or across the 
slope. Define a 1 2 3 system rotated about the 1 axis by 
an arbitrary angle AZ . This is the slope azimuth. The trans- 
formation between these coordinate systems is: 


XI ' 


1 0 0 


XI 

X2 ' 

= 

0 CA SA 


X2 

X3 ' 


0 -SA CA 


X3 


where 


CA = COS (AZ) 

SA = SIN (AZ) 

___ In the primed system, the vehicle has some velocity vector 
V at an angle from the vertical THAVG . Think of this as the 
velocity of the vehicle with respect to the atmosphere. (The 
atmosphere may be moving with some steady wind velocity which 
will be added later.) 

This velocity vector is expressed in the primed system as: 
V(l') = V y = VTERM • COS (THAVG) 

V (2 ) = 0. 

V (3 ) = VHO = -VTERM * SIN (THAVG) 

The atmosphere (or wind) may be moving at any arbitrary 
azimuth direction with respect to these velocities. We assume 
that the wind is blowing horizontally, rather than up or down 
the slope. Thus, we have a vector diagram as shown in figure C3. 


WIND 


3 ” 



Figure C3.- Wind Velocity with Respect 
to Vehicle Velocities in 
the Primed Coordinate System 


(C4) 

(C5) 

(C6) 


(C7) 

(C8) 

(C9) 
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(CIO) 

Ceil) 

(C12) 

These are then the components of the total velocity of the vehi- 
cle with respect to the surface expressed in the 1 ,2 ,3 co- 
ordinate system. 

Now we must derive a relationship between these coordinate 
systems and the vehicle body axes . The vehicle will have a pitch 
angle defined by the initial angle, THAVG, combined suitably 
with a pitch angle produced by a wind gust. Because the wind gust 
is assumed to be in the same direction as VWIND, the vehicle 
will pitch about an axis normal to WIND (namely the p axis 
shown in fig. C3) . To simplify combining the gust pitch angle 
with THAVG, we will assume that both angles can be treated as 
vectors. This assumption is valid for angles up to about 20°. 

THAVG is an angle about the 2' axis. Figure C4 illustrates 
the geometry. 


Vectorially adding the wind velocity to the velocity compo 
nents above, we have 

VOL') = V v 

- V (2 ') = -WIND • SIN(AZWND) 

V (3 ') = VHO + WIND • COS (AZWND) 


AZWND 


3' = Z 



Figure C4 . - Orientation of the Primed Co- 
ordinate System with Respect 
to the Vehicle Body Axes 


Vector addition yields 

TY = THAVG -i- THGST • COS (AZWND) 
TZ = THGST • SIN (AZWND) 

THETA = (TY * TY + TZ • TZ) 


(C13) 

(C14) 

(C15) 
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THETA is the total vehicle pitch angle. Sign is not a prob- 
lem here since AZWND can vary arbitrarily from 0 to 2jt radians. 
THETA will be an angle about an axis rotated AZW from the 2' 
axis, where 


(C16) 
(Cl 7) 

Let us define the 1 ,2 ,3 system with respect to the 

1 , 2,3 system as shown in figure C5: 



SIN(AXW) = SAW = TY/THETA 
COS (AZW) = CAW = TZ/THETA 


Figure C5.- Orientation of the Double- 
Prime Coordinate System 
with Respect to the Prime 
Coordinate System 


Then the transformation becomes 


XI" 

X2 ” ” 



1 

0 

0 

CAW 

0 

SAW 


XI ' 
X2' 

-[cl 

XI 

X2 

X3" 


0 

-SAW 

CAW 


X3' 


X3 


where 



0 

(CA • CAW - SA • SAW) 
(-CA • SAW - CA • SAW) 


0 

(SA • CAW + CA • SAW) 
(CA • CAW - SA • SAW) 


(C18) 


(C19) 
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As stated on the previous page, the vehicle is pitched an 
angle THETA about the l ' ' or P axis. Then, the vehicle 
may be at an arbitrary roll angle PHI about the roll axis. 

The transformation matrix relating the vehicle-body-axis system 
to the system is 


_x b" 


XI ' ' 


CTH 0 -STH 


XI ' ' 

Y b 

- M 

X2 ' ' 

= 

STH • SPH CPH SPH • CTH 


X2 " 

z 


X3" 


STH • CPH -SPH CTH * CPH 


X3" 

b 




- 


*— -* 


where 

CTH = COS (THETA) 

STH = SIN (THETA) 

CPH = COS (PHI) 

SPH = SIN (PHI) 

The beam directions can be expressed as unit vectors in body 
coordinates. For instance, if ALF is the angle about the 

axis to the beam -X^ plane, and if GAM is the beam splay 

angle from the centerline, then the beam direction is 

T = COS (GAM) T v + SIN (GAM) • SIN (ALF) T 
beam X b * b 

- COS (ALF) • SIN (GAM) T 

h 


where 


beam , 



1 , and 1 are unit vectors. 

\ z b 


We can obtain the beam direction in the original inertial co- 
ordinate system (1,2,3) by using the transformations derived 
above. That is, if 

M - M H 


then 


1 


beam 


(1,2,3) 


-M 


T - 


beam 


(W x b) 


where the superscript T denotes matrix transposition. Since 
there are four beams, the Monte Carlo program uses the following 
ope rat i on : 


(C20) 


(C21) 

(C22) 

(C23) 

(C24) 


(C25) 


(C26) 

(C27) 
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!ybm 1 4x3 - [xbm] 4x3 [ai 3x3 


where YBM has rows which are the beam-direction components along 

the 1,2, and 3 axes, respectively, and XBM has rows which 

are the beam-direction components along the X , Y , Z axes, 

. . -i b b b 

respectively. 


Now that we have the beam directions in the 1,2,3 system 
and the normal to the surface slope in the same system, we can 
compute the incidence angle for each beam by: 


TAN (A INC 1) = 


[ YBM(I) 1 X [YSLQP1 
[ YBM(I) ] • [ YSLOP ] 


where X denotes the vector cross-product, the dot denotes the 
vector dot-product, and I denotes the appropriate row of the 
matrix. This method assumes that the same surface slope is seen 
by each beam. 


We must know the velocity along each beam to determine the 
doppler frequency and the receiver sensitivity threshold. 


First, the velocity in the 1^,2^ >3 ^ system is found by 
solving: 


V(1") 

= 

10 0 


v(l') 

V(2") 

= 

0 CAW SAW 


V (2 ) 

V (3 ) 

= 

0 -SAW CAW 


V (3 ) 


Then, 


v v 


V(l") 

*b 



V 

= [B] 

V(2") 

b 





V(3") 

b_ 


_ 


The velocity along each beam component in the system is found 
by combining the dot product of each velocity component in the 
body-axis system with appropriate rows of the XBM matrix. 


The preceeding equations are solved to obtain beam velocities 
and incidence angles for each pass of 1000 passes in the Monte 
Carlo program. For each pass' different values for SLOPE, 


(C28) 


(C29) 


(C30) 


(C31) 
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VWIND, THGST, AZWND, AZ, and PHI are chosen from random- 
number-generator subroutines. The values for THAVG, VTERM, H, 
and radar beam angles are parameters which do not change in a 
given set of 1000 cases. Thus, each pass can be interpreted as 
one possible trajectory point. By making 1000 passes, we obtain 
the statistical information for that trajectory point for 1000 
random geometrical conditions. 


B. Radar Lock Equations 

The Monte Carlo program uses a subroutine that, for a given 
beam velocity, incidence angle, and altitude, will define whether 
the beam is locked or not. Using the beam velocity, the program 
computes the doppler frequency and then computes the receiver 
threshold using curve-fitted equations for the sensitivity curves 
given in reference 4 (for simulation of cw radars). Then using 
the altitude H and the incidence angle, the program computes 
the range to the surface along the beam direction as 

R = H/COS (AINC1) 

This method of computing range is based on the assumption that 
the surface is a plane that slopes at the angle SLOPE used in 
the preceeding geometry equations. This plane is assumed to have 
the same slope at all the points where the beam intersect with 
the surface. Figure C6 below shows the geometry for a single 
direction. Ct is the slope in the direction of the beam direction. 



Figure C6 . - Geometry Used in the Parachute-Descent Static 

Monte Carlo Radar-Lock Program to Compute Range 


(C32) 
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Actually 


R = D/COS (AINC1) 

for this geometry. However, using H instead of D causes less 
than 1.5 dB of conservative error in the radar range equation and 
simplifies the computation for range. The radar range equation 
used is 


PR = C + 4.343 Log e (1/COS 6 * R 2 ) + CSEO 

where P is the power received, CSEO is the sigma zero from 
the reflectivity model, C is a constant depending on trans- 
mitted power, antenna gain, carrier frequency, microwave losses, 
etc, R is the range, and 9 is the incident angle. 

For the ICW radar, an efficiency factor is added to equa- 
tion (C34) , but for a Bessel sideband radar, the range-dependent 
Bessel coefficient must be included. The cross-lobe powers are 
computed as previously described. 

If the power received is less than the threshold, the beam 
is interpreted as being unlocked. If the side-lobe power is 
above the threshold, a potential cross-lobe lock is said to exist. 
The dynamics of the tracker, bandwidth, etc, cannot be included 
in a static point evaluation of this type. 


(C33) 


(C34) 
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INTERLEAVED VERSUS SIDE-BY-SIDE APERTURE ANTENNAS 
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The side-by-side aperture antenna system offers the advantages 
of lower beam coupling, better side-lobe control, greater beam- 
positioning flexibility, greater beam-shaping flexibility, and 
higher redundant reliability. Its disadvantages are increased 
weight and size. 


In the radar performance study, the probabilities of cross - 
lobe and side-lobe lockup were defined in terms of selected re- 
flectivity curves and arbitrarily -assumed side-lobe levels. Be- 
cause increasing the aperture size improves the performance of an 
interleaved antenna system, one of the considerations for all 
modulations was the variation in side-lobe level. The other im- 
portant consideration was the increase in coupling for the re- 
spective channels to the interleaved antenna system. For all 
cases considered, we assumed that the planar array systems would 
be comprised of rectangular slots which are not space -tapered . 


In most analyses of the coupling of the rectangular slot, 
Babinet's theorem is invoked to prescribe the analysis in terms 
of a dipole antenna in free space (no ground plane). Using the 
electric and magnetic vector potentials, the electric and mag- 
netic fields are given by 


E 


i W ^ ] 

- jwA - 1 - V (v * A) - -VXF 



k 2 




H 


- jwF - -l— y(y • F) + — VXA 

* *“ 


K 2 M -1 


(Dl) 


For the rectangular slot shown in figure Dl, the electric 
field is given by 


E - --VXF = - 


er sin 


a r ra 6 r ( sin e ) a r 


S_ <L 

Sr Se 


F rF Q 
r 9 


0 

0 


(D2) 
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or 


: - £ [u* e y 


where 


F cos 6 
z 


F sin 0 
z 


/ vmL exp \ 

\ 4itr ) 

_ / -VmL exp" jkr \ 
\ 4«r ) 


(cos 6) 
(sin 0) 


or 


- ( — g~ ikc ) (y + y 


and for 


H = -jwF - 


V^V • 


where in spherical coordinates 


V * F = 


F s in 6 

u 


) 


■ ^[s^ 4 F r )j + 

( v ' 0 - h[y ■ ; ) * r _ + r j h [(* - F ) *6 } 


then 


H 

r 


(VmL 

-jkr\ 

exp 

(cos 8) j 

f 2 | 2 \ 


V 

2 it ) 

^Z Q r 2 jwpr 3 / 



(VmL 

ex P~^ kr ^ 

(sin 0) j 

(jJs_ + _L_ + 1 \ 


V 

) 

y 0 z Q r j wpr y 


(D3) 


(D4) 


(D5) 


(D6) 


(D7) 

(D8) 


(D9) 
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Utilizing Babinet's principle. 


~ * 

H = E 
o o 


(DIO) 


E = -Z H 
o o o 


(Dll) 


Then, the fields of the slot ( E cp s » H g s ^ can be re P resented b Y 
the fields of the wire as 


E cps H 9 w 


6s z 2 

o 


(D12) 


For the two dipoles shown in figure D2 , the basic array expres- 
sion can be utilized: 


E „ < e ^> ■ r*r 


) (*«>,¥)) E A n * [ E(e ’ W ] \ft <*>] <D13) 


where 


. N-l 

?<«>-E 


A exp~ 
n 


and for a uniform array 


(N-l)/2 


E 


A n cos (nv + 9n) 


or for the 2 -dipole array 


E ( 0, v) = 2 Ed (O', • ) 


r kd cos 0 - \|r 
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where 


Ed(0,Cp) = Id 


cos 


/ cos 6 \ 

\ *—r) 


sin Q 


(D17) 


The coupling is normally defined in terms of the mutual impe ■ 
dance. Considering identical antennas |Zn = Z 22 j 


"Vi' 


Z11 

2 12 


M 

V 2 _ 



■ i 

C\l 

1 — l 
N 


1 

M 

ro 


and for 


we have 


Z 12 = jZ 12 | 


(D18) 


(D19) 


W, 


W 2 


ill 

[Rix + 

1 Z 12 1 

COS 4* 

ill 

[r 22 + 

|Zi 2 | 

cos fa - 


(D20) 


Then, the total power W of the antennas is 

W = Wi + W 2 = 2 | Ix| 2 ^Rii + (R 12 ) cos 


(D21) 


and 

E(e,cp) = 



/jt cos 6\ 


• « 

2 I « 1* 

\ 2 ) 

cos' / 


/kd cos cp-Tlr v 

| 2 [R 11 + (R 12 ) ^cos j 

sin 6 


[ cos \ 2 ) 


D22) 


and the gain of the two dipoles compared to a single dipole is 


G 

a 

E 2 
a 

2(Rn) 

/ kd cos Cp - \lr ' 

[cos- ( a 

G d = 

E d 

R 11 

+ (R 12 ) /cos \|r_^ 


(D23) 
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Now, since the power transmitted is 

W t = if [Re(Z u )] 

and the received power in a matched load is 


J r = if [Re(Z 22 )J = 


[Re(Z 22 )] 


then the field coupling C is 


C = jl = I 

W t 4 (Rii)(R 22 ) 


where 


Rn 


Re(Zn) 

R 2 o 


Re (Z 22 ) 


Because the coupling may be reduced by orientation of the 
inclination angle (see fig. D3), the coupling C can be expressed 
(see ref. 15) in terms of the directivities D^_ ) as 

C - [ C „( R o)] [ D t (\)] [“rCr)] ( N qt , qr ) < D; 

for a rectangular aperture with dimensions (a,b) 
where 


D ( C P) = (cos" q>. 


/cos uV 


U - IT 


(t ) sl1 


(D'iO) 
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V = « cos cp (D31) 

F ~ 1 when ^ and t < 1 (D32) 

Curves illustrating the relative magnitude of the coupling co- 
efficient c ( r q ) and ttie directional constants ^ and D 

(9.) are shown in figures D4 and D5. Figure D4 shows that C ( R Q ) 

decreases at the rate of 6dB/octave # In figure D5, Curve I shows 

the case where both <P and 9 =0°, and Curve II shows the 

t r 

case where both 9 fc and 9 r =Jt/2. As shown in figure D4, for 
case I ^ when 9^ and 9 fc = 0°^, the isolation between inter- 
leaved antennas would be about 16 dB for R Q /^ = 0.5 and 22dB 
for R^y/A = 1*0. For the side-by-side aperture antennas, the 

curves can be extrapolated (using 6dB/octave) to greater than 60dB 
for separate 18x18 -in. apertures, if we make the assumption that 
the effective array R^ >18 in. 

From equation (D13), we can see that if the number of ele- 
ments N of the interleaved array equals the number of elements 
of an equivalent array without interleaved elements, then the 
array field E^(6,9) is unaltered if E(6,9) is unaltered. 

From equation (D22), the effects of coupling or mutual impedance 
are shown to decrease the gain of the individual element E ( 6 , 9 ) 
by the R 12 term. From these equations, it would appear that 
the basic side-lobe structure is not altered with respect to the 
main beam if the above conditions are retained. The additional 
elements of the interleaved array appear as an equivalent yagi 
antenna (fig. D6); and the excitation in the reference array 
elements produces a traveling wave in the corresponding elements 
of the interleaved array. The current driving the nth element 
is 

I ~ ik x ns \ 

I n <Z') = I <z') \ exp / (D33) 

and the magnetic vector potential in the far field is 
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Figure D5 . - Broadside Coupling of Open-Ended Slots as a Function of the Orientation 
of Two Slots in a Planar -Array Antenna 
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If we recognize the effect of the interleaved array to be the 
ex P (~ J k x s where 


cos 


k^s = (cos ks) - ^exp 5 a \ 


(D35) 


and 


(h,a) = 


(cos kn) (saJsL) - a (cos kh) 

sin kh - kh (cos kh) 


(D36) 


then we can see that the change in the far field pattern depends 
on: the equivalent slot dimensions (a,b) related to the dipole 

dimensions (h,a); the slot spacing (ks) or phase relation- 
ship; and the excited element with respect to its location on 
the array. If we note that the exponential term, exp^-jk s^, 

can be expressed in an alternating power series that can be maxi- 
mized by selecting the spacing S in multiples of A/2, then 
we see that the maximum contribution from the other elements would 
occur at the center of the array, with a taper at the ends of the 
array. For a mul tiple -beam array, the main beams would not occur 
at the center, and the side -lobe structure would be increased on 
the sides of the main beams near the center of the array. 
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MARTIAN SURFACE REFLECTIVITY 
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Two sets of reflectivity curves were used in the simulation. 
One set consisted of two upper and lower bound curves with the 
corresponding mean curve as described in reference 8. The other 
set of reflectivity curves that was used is based on the following 
assumptions: a certain portion of the local terrain will have a 

specular return similar to that measured on earth (ref 16) ; the 
greater portion of the return will be from the diffuse surface, 
which means that this return will be dependent on the local ter- 
rain roughness and electrical characteristics; and the beamwidth- 
illuminated area will tend to average these variations into an 
effective mean curve at higher altitudes. 


For the diffuse surface components, the normalized radar cross- 
section per unit area oo can be expressed in terms of the 
Fresnel reflection coefficients. 


r + = (Y g cos 6) - (Y 2 - sin 2 6)^ 
(Y 2 cos 6 ) + (Y 2 - sin 2 0) ^ 


(El) 


where 


Y = 


(W v) ! ” ( e / £ °) 


% 


(E2) 


Then, if a curve-fitting model such as Beckmann’s model is used 
for the microscopic roughness from reference 17 (fig. El) 


P(6) 


1 

(cos 4 6 + R sin 4 6)^^ 



erfc (K cot 0) 


•>]i) 


(E3) 


The effect of surface slope and beamwidth averaging can be 
incorporated into the sigma zero expression 

( c o) K ■ c 0 { e *e h" 2 <p o 1 ( s ! 1 s i)fe)]| 

where it can be seen from figure E2 that 

2 h 

tan cp = — cos 0=2 tan a cos G 


(E4) 


(E5) 


where tan a is the surface height divided by the surface decor- 
relation distance, € is the incidence angle to a flat surface, 
and the integration intervals are 


245 



APPENDIX E 



Incidence angle, deg 


Figure El.- Beckmann Reflectivity Curves 
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Figure E2.- Effect of Surface Slope and Beamwidth Averaging 
on Mars Surface Reflectivity 


a = H 
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(E7) 


If we consider the microscopic roughness (a < X) to be super- 
imposed on a microscopic rough surface (a > A), the effects of 
the beamwidth smoothing can be seen from figure E3 . In calculat- 
ing the effects of beamwidth averaging, the surface roughness was 
considered to be normally distributed 


R h (M) = 



(E8) 


with slopes 


R tan 9 <5,,1> 


tan 2 cp 


exp 



(E9) 
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Figure E3 . - Mean and Variance Curves for Mars Surface 
Reflectivity at an Altitude of 4000 ft 
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p(tan cpi,tan cp 2 ) 



' 1 



(tan 2 cpj |l - [p(S,n)] 2 |* 

cpi + tan 2 cp P - 2 (tan cp n ) ( tan cp 0 ) [ 0 ( £ 
2 tan 2 cp o |1 - [p(|,T])] 2 | 



Now we can define the exponential autocovariance function 


= tan 2 cp Q 


exp 



and the sample average 
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r b 
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Now we evaluate the variance terms ( s f > s f) of the imaginary and 

real components (vertical and horizontal) of the surface slopes 
at CP]_ and cp 2 < rt/4 radians. 
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RANGE SIGNAL BANDWIDTH 
FOR THE MODIFIED LM RADAR SYSTEM 
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The geometry used in the derivation of the range signal band- 
width is shown in figure FI. 



Figure FI.- Frequencies at the 3-dB Points of the 
Range -Antenna Beam in the i - ~E 
Plan Z 


The figure indicates the frequencies at the 3-dB points of 
the range antenna beam in the plane defined by the vectors 
i^ and B. These frequencies are proportional to the range and 

velocity of the beam. Equation (FI) is used to compute the band- 
width in this plane: 


EW 


f RZ + f D2 


f Rl + f Dl 


(FI) 


where BW is the range signal bandwidth, f R is the frequency- 
proportional to the range, f R is the frequency proportional to 
the velocity, and where f R2 is always greater than f R ^ . 


Note that is the Doppler BW in the i - F plane is 

u 


bw d = 


D2 


f 


Dl 


and if the BW proportional to the range is 


(F2) 
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BW r 


R2 


f 


Rl 


then equation (FI) can be written as 

BW = BW d ± BW 
R D 


The Doppler BW will be additive when the velocity vector is 
above the conic generated by swinging the beam-pointing vector 
about the vertical, and subtractive when the velocity vector is 
within the cone. 


(F3) 


(F4) 



Figure F2 . - Geometry for Calculating the Doppler 
Bandwidth BW 

K 


From figure F2 we can see that d is proportional to BW R . 
d can be found from 




h_ 

cos ^cp + 




Then , 


BW d » (SF) R0 tan cp 
R 


(F5) 


(Fb) 


where SF, the range scale factor, is in units of Hz/ft. 

Figure F3 shows the geometry used to obtain the Doppler band- 
width, BWr, in the plane described in figure FI. 
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Figure F3.- Geometry for Calculating the Doppler 
Bandwidth BW^ 

In this figure, is the range-beam pointing vector, V. 

is the velocity vector, and i , i , i are the unit vectors 
along the inertial axes. 1 L 

The problem is to determine the component of V. perpendicular 

to B^ and the component in the plane defined by T and B\ 

Z i 

Since the range beam points along the U body axis 

B. = [A] T [B] = 

where A is the inertial-to body-axis transformation and B is 
the range-beam pointing vector in the body -coordinate system. 

A unit vector perpendicular to T and B. is given by 

Li X 



0 



L x 


"- A 12 

X ' 

1 

p = r z x b i = 

0 

X 

A 12 

r x 

= 

A 11 

h 


1 

1 N 

1 


A !3 



0 

0 


We want to find the unit velocity perpendicular to B. in 
the plane defined by T and B.. This is 1 

u 1 


hhx 

A 12 L Y 

A 13 L Z 


(F7) 
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where y is the angle between the range -beam pointing vector 
and the velocity vector. Otherwise, the bandwidth defined by 
equation (FI 3 ) is used. 
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